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HELP,  A  MULTI -MATERIAL  EULERIAN  PROGRAM  FOR 
COMPRESSIBLE  FLUID  AND  ELASTIC-PLASTIC 
FLOWS  IN  TWO  SPACE  DIMENSIONS  AND  TIME 


ABSTRACT 

Volume  I  details  a  new  numerical  model  for  solving  multi-material 
flows  which  are  functions  of  two  space  dimensions  and  time.  The  program, 
called  HELP,  is  basically  Eulerian  but  also  contains  Lagrangian  features 
for  the  explicit  definition  of  interface  positions.  The  program  is  gen¬ 
eral  in  the  sense  that  any  number  of  materials  can  be  present  in  a  given 
Eulerian  cell  and  no  special  difficulties  arise  for  flows  involving 
extreme  material  distortions. 

The  program  is  suitable  for  compressible  media  with  strength  depen¬ 
dence.  Compressibility  is  included  by  an  equation  of  state  giving  pres¬ 
sure  as  a  function  of  density  and  specific  internal  energy.  Strength  is 
included  by  means  of  a  constitutive  equation  giving  the  deviatoric 
stresses  as  fimctions  of  the  elastic  and  plastic  strains. 

In  this  Volume  I  report,  the  basic  equations  are  developed,  the 
computer  program  is  documented  and  sample  results  are  given  from  six 
applications . 
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I.  INTRODUCTION 

HELP  is  an  attempt  to  address  and  resolve  certain 
familiar  and  long-standing  problems  in  hydrodynamic  computing 
which  are  not  adequately  resolved  by  pure  Lagrangian  or 
Eulerian  methods.  Describing  the  motion  of  a  material  inter¬ 
face  in  a  Lagrangian  code  is  automatic,  but  the  numerical 
solution  becomes  unsatisfactory  if  the  material  undergoes 
extreme  distortion.  In  Eulerian  codes,  on  the  other  hand, 
extreme  distortion  poses  no  problem,  but  the  accurate  treat¬ 
ment  of  surface  motion  (i.e.,  a  material  interface  or  a  free 
surface)  requires  a  detailed  knowledge  of  surface  location 
which  is  not  contained  in  a  purely  Eulerian  program.  .  In  the 
present  work,  massless  tracer  particles  are  introduced  which 
define  the  surface  positions  and  move  across  the  Eulerian 
grid,  so  that  the  program  yields  a  Lagrangian- type  definition 
of  the  moving  surfaces.  At  the  same  time  the  code  does  not 
lose  the  Eulerian  capability  for  treating  extreme  distortions 
without  difficulty.  Transport  from  Eulerian  cells  containing 
either  free  surfaces  or  interfaces  between  two  or  more 
materials  is  done  in  a  manner  which  takes  account  of  these 
surface  positions. 

The  HELP  code  has  evolved  from  four  major  hydrodynamic 

programs  developed  over  the  past  fifteen  years.  These  codes 

have  been  used  in  a  wide  range  of  applied  problems  in  shock 

hydrodynamics.  Starting  from  the  one-material  code  RPM^^^ 

the  HELP  development  reported  here  adds  the  capabilities  of 

describing  multi-material  interactions  and  treats  material 

strength  as  elastic-plastic  rather  than  rigid-plastic.  RPM, 

r  2 1 

in  turn,  was  an  extension  of  the  OIL  compressible  fluid  code'-  ■' 
The  name  HELP  means  Hydrodynamic  plus  Elastic  Plastic. 
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to  include  material  strength  in  the  r igid-per£ectly  plastic 

r31 

approximation.  OIL  in  turn  evolved  from  the  PIC  ^  ■'  code  and 
was  developed  by  replacing  the  discrete  PIC  particles  by  a 
continuum.  This  latter  transition  to  a  continuum  calculation 
was  made  in  order  to  permit  the  treatment  of  very  small  com¬ 
pressions  (such  as  those  usually  accompanying  processes  in 
which  material  strength  is  important),  without  having  to  use 
an  extremely  large  number  of  particles  for  adequate  mass 
resolution.  A  disadvantage  of  this  transition  was  the- loss 
of  certain  Lagrangian- type  features  of  PIC,  and  specifically 
the  ability  to  treat  flows  containing  more  than  one  material. 
The  major  purpose  of  the  work  on  HELP  (other  than  the  effort 
to  treat  strength  with  elastic  and  plastic  effects)  has  been 
the  restoration  of  this  multi-material  capability,  while 
retaining  the  ease  of  treating  small  compressions  and  the 
computing  economy  of  a  continuous  Eulerian  model. 

In  all  of  the  programs  mentioned  above,  the  basic 
difference  equations  are  explicit  and  are  accurate  to  terms 
of  first  order.  The  difference  equations  used  in  HELP  will 
be  given  in  this  report,  but  the  interested  reader  is 
referred  also  to  Refs.  1  through  3,  and  literature  cited 
therein,  for  further  discussions  of  the  properties  of  these 
or  very  similar  difference  equations. 

The  present  code  development  has  successfully  led  to 
the  ability  to  treat  extreme  strength-depende'nt  deformations 
in  multi-material  interactions.  Examples  are  the  axisymmetric 
interaction  of  projectiles  of  various  shapes  with  target  plates 

[41 

of  different  materials. ■■  Similarly,  in  Section  VI  of  this 
report  are  seen  results  of  calculations  from  a  half-dozen  re¬ 
presentative  problems.  It  is  still  possible  with  the  program 
to  treat  the  simpler  cases  of  the  earlier  codes;  specifically, 
strength  effects  can  be  ignored  by  simply  bypassing  one  phase 
of  the  calculation,  and  no  changes  are  required  to  apply  the 
program  to  flow  fields  containing  a  single  material. 
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HELP  cannot  be  utilized  effectively  without  investing 
a  considerable  effort,  perhaps  1  to  2  manpower  for  several 
months,  in  becoming  familiar  with  the  basic  numerical  method, 
the  details  of  the  FORTRAN  program,  and  in  gaining  computing 
experience  with  the  code.  While  every  effort  has  been  made 
to  make  the  code  as  reliable  and  general  as  possible,  it 
should  be  noted  that  any  new  application  may  give  rise  to 
problems  or  anomolous  results  (e.g.,  a  failure  of  the  pressure 
iteration  in  a  mixed  cell)  which  can  effectively  be  diagnosed 
and  resolved  only  by  persons  who  are  sufficiently  familiar 
with  the  entire  program  that  they  can  rapidly  identify  and 
remedy  difficulties  with  a  calculation.  Some  such  general 
statements  are,  of  course,  true  of  any  2-D  hydrodynamics  code 
as  applied  to  a  new  problem.  In  the  case  of  HELP,  working 
familiarity  with  one  of  the  predecessor  codes  RPM,  OIL  or  PIC 
will  be  of  considerable  value  in  reducing  the  above  learning 
period.  In  the  present  report,  an  attempt  is  made  to  give 
the  potential  user  a  comprehensive  introduction  to  both  the 
numerical  methods  and  the  computer  program. 

Some  comments  will  guide  the  reader  through  the  text 
of  this  report.  The  next  section  (II)  is  a  general  descrip¬ 
tion  of  the  numerical  method.  In  Section  III,  the  treatment 
of  strength  with  the  elastic-plastic  constitutive  equation, 
is  reported  in  some  detail  because  these  results  are  essen¬ 
tially  new.  Similarly,  Section  IV  describes  in  some  detail 
the  method  which  is  used  for  propagating  interfaces  and  free 
surfaces,  and  the  manner  in  which  information  on  interface 
position  is  used  in  affecting  the  Eulerian  hydrodynamic  solution. 
Section  V  details  the  pressure  iteration  which  is  necessary 
in  order  to  determine  the  pressure  (and  material  densities) 
in  cells  containing  more  than  one  material.  Section  VI  pre¬ 
sents  sample  results  obtained  with  HELP.  Section  VII  is  a 
guide  to  the  FORTRAN  program,  giving  such  information  as  the 
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flow  diagram,  the  functions  of  the  various  sub^routines ,  the 
input  for  generating  and  restarting  problems,  and  the  dictionary 
of  FORTRAN  variables.  Volume  II  contains  a  listing  of  the  FORTRAN 
program.  The  program  is  interspersed  with  comment  cards, in  an 
attempt  to  facilitate  learning. 

Learners,  users  and  improvers  of  HELP  are  invited  to 
communicate  with  the  authors  where  such  improvement  might, 
facilitate  their  efforts  or  lead  to  improvements  in  the 
program. 
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II.  GENERAL  DESCRIPTION  OF  THE  NUMERICAL  METHOD 

In  the  present  section  we  give  the  conservation  equa¬ 
tions  and  the  numerical  model  which  is  used  to  treat  these 
equations.  The  equation  o£  state  which  is  used  for  condensed 
media  is  also  presented  with  equation-of-state  parameters  for 
some  nineteen  materials. 


2.1  CONSERVATION  EQUATIONS 

Space  is  divided  into  fixed  Eulerian  cells  through 
which  the  fluid  moves.  To  arrive  at  expressions  for  the 
rate  of  change  of  total  mass,  momentum  and  energy  within 
such  a  cell,  it  is  convenient  to  start  with  the  equations 
of  motion  in  the  form; 


ie.  = 

at 


(1) 


p 


(2) 


C3) 


Here  is  the  stress  tensor,  which  can  be  regarded  as 

the  sum  of  the  hydrostatic  stress,  "S^jP,  and  a  stress 
deviator  tensor,  ,  i.e.. 


(4) 


and  E^  is  the  total  energy  (kinetic  plus  internal)  per 
unit  mass.  Tensor  notation  is  implied,  so  that  repeated 
indices  denote  summations. 
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Expanding  the  convective  derivatives  in  Eqs.  (2)  and 
(3),  D£/Dt  =  3f/3t  +  3£/3x^,  then  adding  Eq.  (1)  times 

Uj  to  Eq.  (2),  and  Eq.  (1)  times  E,j,  to  Eq,.  (3),  and 
collecting  terms,  gives 


k  (o^j)  '  ^  “ij  -  air  (55 

k  W)  ■  kr  -  air  (p“A)- 


For  the  developments  to  £ollow  it  is  Idesirable  to 
replace  these  di££erential  equations  by  the  analogous  integral 
equations,  obtained  by  . integrating  over  the  cell  volume,  V, 
and  then  converting  the  volume  integral  o£  divergences  to 
sur£ace  integrals  over  the  cell  sur£aces.  Equations  (1), 

(5) ,  and  (6)  then  become 


J  pdV  =  -  jT  pu^n^ds 

V  s 


(7) 


Ft/  PUjdV  =  /ajjn^ds  -  /  pu.u.n^ds 


(8) 


/  pE^dV  =  /o..u.njds  -  J  pUjE^njds  . 


C9) 
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2,2  DIFFERENCE  EQUATIONS 


2.2.1  Division  into  Phases 

It  is  convenient  to  express  the  integral  conservation 
relations,  Eqs,  (7)  through  (9),  as  finite  difference  equa¬ 
tions  over  the  time  step  At  and  also  to  decompose  the  total 
stress  a^j  into  its  deviator  and  hydrostatic  components, 
according  to  Eq.  (4).  This  gives,  for  the  increments  of 
total  mass  (m) ,  momenta  (mu.j  and  energy  (mE^)  within  the 
cell 


Am  = 


f 


-At  I  pu^n-ds 

(10) 


(mu.)  =  -  At  i 

r 

Pn.ds  + 

'  r  J 

J 

J  s^jn^ds  -  At/*  (pu^UjJn^ds 

iff*  f* 


A^mE^j  =  -At  ^  Pu^n^ds  +  At  J  s^^u^n^ds  -  At  J"  ^pu^E,p)n^ 


(11) 


ds 


(12) 


Here,  the  terms  on  the  right  are  divided  into  increments  due 
to  pressure  forces  on  the  cell  surface  (first  column),  those 
due  to  the  stress  deviator  forces  on  the  cell  surface  (second 
column)  and  the  increments  (third  column)  due  to  the  trans¬ 
ports  of  mass,  momentum  or  energy  through  the  surface  of  the 
cell.  These  three  types  of  contributions  are  accounted  for 
in  the  computation  in  distinct  phases.  Specifically,  during 
each  time  step  all  cells  are  updated  for: 


Pressure  effects  in  Phase  1 

Effects  of  the  stress  deviators  in  Phase  3 

Transport  effects  in  Phase  2 
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(The  calculations  are  performed  in  the  order  listed,  naming 
of  the  phases  having  to  do  with  the  history  of  their  develop¬ 
ment.)  In  the  discussion  that  follows,  the  calculation  of 
the  terms  on  the  right  of  Eqs.  (10)  through  (12)  are  described 
sequentially  starting  with  Phase  1. 

Some  preliminary  definitions  will  be  useful.  Super¬ 
script  n  on  a  variable  refers  to  the  value  of  the  variable 
at  the  beginning  of  the  time  step  and  superscript  (n+1)  denotes 
the  value  at  the  end.  In  this  discussion  a  typical  cell  in 
the  interior  of  the  grid  is  considered;  Section  2.4.1  has  a 
discussion  of  the  special  conditions  which  describe  the  cal¬ 
culation  at  the  grid  boundaries  or  at  the  axis  of  symmetry. 

For  a  typical  cell,  denoted  by  a  value  of  the  index  k,  the 
dependent  variables  for  that  cell  are  written  P(k),  u(k)  , 
v(k),  Ej^(k),  m(k),  representing  respectively  ;the  pressure, 
radial  and  axial  components  of  velocity,  the  'specific  internal 
energy  and  the  mass  for  cell  k.  The  adjacent  cells  above, 
below,  to  the  right  and  left  of  k  will  be  designated 
respectively  as  ka,  kb,  kr,  and  kl .  Here  the  terms  above, 
below,  right  and  left  refer  to  a  cross-section  view  of  the 
cells,  in  which  the  left  border  is  parallel  to  the  axis  of 
symmetry  with  z  increasing  upward  (see  Fig.  1).  Each  cell 
is,  then,  the  torus  obtained  by  rotating  the  'rectangle  (since 
Ar  ^  Az  in  general)  about  the  axis  of  symmetry. 

2.2.2  Phase  1 — The  Effects  of  Pressure 

In  the  following  discussion,  the  Phase  1  calculation 
is  described  first  for  one-material  cells.  The  changes  made 
to  extend  the  Phase  1  calculation  to  a  multi-material  cell 
are  then  given  in  a  concluding  paragraph. 
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2 . 2 . 2 . 1  Continuity  Equation,  (10) 

No  contribution  to  Phase  1. 

2. 2. 2. 2  Equation  of  Motion,  (11) 

1.  Axial  Motion:  In  this  case,  referring  to  the 
pressure  integral  in  Eq.  (11),  u^  =  v  and  pj  =  -1,  0,  +1 
are  the  axial  components  of  the  unit  normal  to  the  bottom, 
sides  and  top  of  the  cell  respectively.  Equation  (11) 
therefore  gives  for  the  Phase  1  contribution 

A^(mv)  =  P^TT^r^  -  r^^At  -  P^^TT^r^  -  r^jAt 

where  r  and  r  are  the  radii  of  the  outer  and  inner  cell 
2  1 

surfaces,  respectively,  and  ,  P^  are  the,  pressures  on  the 
bottom  and  top  cell  faces.  These  pressures  iare  obtained  from 
the  initial  (time  n)  values  of  cell  pressures  by  a  simple 
average  of  the  pressures  in  the  adjacent  cells: 

P  =  P^(k)  -H  P^(kb) 

^b  2 

P  -  P^(k)  +  P^(ka) 

^a  2  • 

Making  these  substitutions  gives  the  Phase  1  increment  of 
axial  momentum 

2.  Radial  Motion:  To  arrive  at  the  padial  equation 
of  motion,  it  is  useful  to  consider  a  volume  with  the  full 
cell  dimensions  Az  and  Ar  in  the  z  and  r  direction, 
but  with  a  small  angular  dimension  A9  instead  of  2tt  (see 
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Fig.  2).  Then  one  can  compute  the  radial  motion  from  Eq.  (11) 
since  the  motion  corresponds  to  a  fixed  direction  in  space. 

For  use  in  Eq.  (11),  then,  u^  =  u,  and  n^  =  -1,  1,  -A0/2,-Ae/2, 
0,  0  correspond  to  the  radial  components  of  the  unit  normal 
to  the  inner,  and  outer  surfaces,  the  two  side  surfaces  and 
the  top  and  bottom,  respectively.  The  pressure-integral 
contribution  in  Eq.  (11)  becomes 

A  A 

■s—  A  (mu)  =  P„r  A0  Az.  At  -  P^r  A0  Az  At  +  P  Ar  Az  A0  At 
27r  1  ^  1  r  2  s 

where  P^^,  P^  and  P^  are  the  left,  right  and  side  face 
pressures  respectively  and  are  defined  in  terms  of  time  n 
cell  pressures  as  follows. 

.  P"(k)  ♦ 

2 

_  P^(k)  ^  P^(kr) 
r  2 

and 


+ 


p"CkJl)  +  2P^(k) 

r  - 


In  terms  of  these  cell  pressures,  and  using  Ar  =  r^  -  r^, 
the  equation  for  radial  motion  becomes 

AZ  At  . 
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Fig.  2- -Element  o£  volume  for  discussion 
of  radial  motion. 
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2. 2. 2.3  Energy  Equation  (12) 

Computing  the  Phase  1  energy  equation  departs  somewhat 
from  the  standard  procedure  which  is  used  to  compute  the  other 
integrals  on  the  right  in  Eqs.  (1)  through  (12).  Specifically, 
rather  than  compute  the  gain  in  total  cell  energy, 

■"K  *  •  by  integrating  the  pressure  forces  on 

the  cell  surface,  the  original  equation  is  reformulated  to 
give  a  volume  integral  expression  for  the  gain  in  internal 
energy,  the  kinetic  energy  increment  having  already  been 
fixed  by  the  preceding  calculation  of  the  velocity  increments. 
The  resulting  equation  is  preferred  because  it  is  a  simpler 
expression  for  E^^  and  because  it  has  been  used  in  most  of 
our  calculations  and  has  led  to  satisfactory  results.  Starting 
with  the  Phase  1  part  of  Eq.  (12), 


Pu^n^ds  . 


and  since  the  cell  mass  m  is  constant  in  Phase  1  and 

Et,  =  E^  +  i  u-u.  one  can  write 
T  m  2  1  1 


A  (mE  1  +  u-A  (mu.\  =  -  /  Pu-n-ds 

i'  ml.  iiV  1/  J  11 


•'v  " 
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where  P  and  in  the  last  line  are  average  values  over 

the  cell  volume.  But  from  the  Phase  1  momentum  equation 

\(“i) 

or 

SO  that  the  last  terms  on  each  side  of  the  energy  equation 
cancel  and  the  resulting  expression  for  the. gain  in  internal 
energy  is 


where  is  obtained  by  averaging  the  time  (n) 

(pre-Phase  1)  and  time  (n+l)  values  obtained  from  the  Phase  1 
momentum  calculation,  i.e., 
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and,  for  example,  the  time  (n)  value  of  this  divergence  is 
determined  from  the  calculation 


3u . 
1 

1 


3 


(ru) 


+ 


3v 

3z 


,  I  .  v"fkal  -  v"fkb) 

3x^  I  r^  y  iAr  )  2 A z 


In  which  r  =  (r  +  r  )/2,  r  =  r  +  Ar/2,  r,  =  r  -  Ar/2 
c  \  1  2'  r  2  ^  y 

the  radii  of  cell  centers  for  the  central,  right  and  left 
cells  respectively. 


are 


This  completes  the  discussion  of  the  Phase  1  calcula¬ 
tion,  for  a  typical  cell,  of  the  momentum  and  internal  energy 
increments  due  to  the  pressure  forces  on  the  cell.  The  two 
momentum  increments  are  used  to  determine  the  velocity  incre¬ 
ments  by  dividing  by  cell  mass,  and  one  has,  at  the  end  of 
Phase  1,  new  values  of  Ej^(k),  uCk),  v(k). 


2. 2. 2. 4  Extension  of  Phase  1  to  Mixed  Cells 


The  changes  required  to  extend  Phase  1  to  multi-material 
cells  are  minor.  The  velocity  increments  are  determined  as 
before,  it  being  assumed  that  all  the  materials  in  a  mixed 
cell  have  a  common  velocity  (u,  v) .  The  total  increase  in  the 
internal  energy  within  the  cell  is  computed  as  before,  and 
this  energy  is  partitioned  among  the  materials  by  assuming 
that  energy  increments  are  proportional  to  fractional  volumes 
occupied  by  the  various  materials.  These  fractional  volumes 
are  computed  by  dividing  the  mass  of  each  material  by  its 
density.  Densities  and  masses  are  updated  in  the  equation 
of  state  and  Phase  2  routines,  respectively. 
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2.2.3  Phase  5 — ^The  Effects  of  Deviatoric  Stresses 

2 . 2 . 3 . 1  Continuity  Equation,  (10) 

No  contribution  in  Phase  3.  j 

2. 2. 3. 2  Equation  of  Motion,  (11) 

1.  Axial  Motion;  The  Phase  3  contribution  from 
Eq.  (11)  is 

A  ( mu . 

3  '  J 

and  for  the  axial  motion 


s . -n- ds 
ij  1 


Here  s^jn^  is  the  axial  component  of  the  total  stress  on 
a  surface.  For  a  torus  with  rectangular  section,  s^^n^  on 
the  various  faces. is 

^zz 

-  s^  cell  bottom 

zz 

outer  cell  surface 
r  z 

1 

a 

-  s  inner  cell  surface 

r  z 

where  zz  and  rz  subscripts  denote  normal  and  shear  stresses. 
The  top,  bottom,  right  and  left  surfaces  are  indicated 
by  t,  b,  r,  and  £  respectively.  The  equation  for  axial 
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motion  is  therefore 

A^Cmv)  =  -  s^j7r(r^  -  rj  )At  +  Zirr^  Az  At 

-  2'fTr  Az  At  . 

rz  1 

2.  Radial  Motion:  For  the  radial  motion,  consider 
the  element  of  mass  depicted  in  Fig.  2.  Then  in  the  Phase  3 
part  of  Eq,  (11) 


and 


( mu . ) 

=  At  i 

J 

s . - n . ds 
ij  1 


u  • 


u  . 


The  radial  components  stresses  on  the  various 

faces  of  the  mass  element  are 


s 


t 


-  s 


s 


s 


-  s 


-  s 


zr 

b 

zr 

r 

rr 

I 

rr 

00 

00 


top  of  mass 
bottom 
right 
left 

(%A0)  front 
(%.A0)  back  . 


element 
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The  mass  of  the  element  is  mAB/Zir  where  m'  is  the  cell 
mass.  Then  , 


Ae 

2tt 


A  (mu) 

3 


b  \ 

/  2 

2  \ 

1  s  „ 

■  s  ) 

-  r  1 

\  zr 

zr  / 

\  2 

1  / 

Ae 

T 


At 


S  T 
rr  :  2 


A9  Az  At 


-  s^^  r^A9  Az  At  -  SggAz  Ar  Ae  At 
or  (multiplying  by  27r/Ae)  the  equation  for  radial  motion  is 
A^Cmu)  =  Ztt  -  S^^)(r“  -  r!" ) 

+  sl^  -  Sg^Az  Ar]  . 

2. 2. 3. 3  Energy  Equation,  (12)  : 

In  the  Phase  3  term  from  Eq.  (12), 

AglmE^)  =  Aty  s-^.Ujnj_ds  ,  ' 


Si^UjUf  is  the  work  rate  per  unit  surface 
the  various  faces  of  the  torus,  is 


area,  which  for 


(s 
-  (s 


t 

zz 

b 

zz 

r 

r  z 

£ 

rz 


V 


V 


V 


V 


t 

b 

r 

£ 


uM 

zr  / 

cell  top 

uM 

cell  bottom 

zr  1 

^  u^l 

rr  1 

outer  cell 

surface 

£  £x 

rr“ 

inner  cell 

surface 

where  t,  b,  r,  £  denote  stresses  and  velocities  at  the  top, 
bottom,  right  and  left  cell  faces  respectively.  The  interface 
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velocities  are  determined  from  a  simple  average  of  the  time 
n  velocities  in  the  two  cells,  e.g., 

„t  _  v'^Ck)  +  v'^Cka:) 

V - ^  2 


£  _  u^(k)  u^Ckl). 
u  2 

The  equation  for  the  Phase  3  change  in  total  cell  energy  is 

A  (mE„  ]  =  (s^  v^  +  -  s^  v^  -  s^  u^W(r^  -  r^\At 

3 \  T  /  \  zz  zr  zz  zr  /  V  2  i  / 

+  (s^  v^+  s^  u^]  2it  r.  Az  At 
\  rz  rr  /  2 

I  I  I  ^  I  l\o  A 

-  (s^,v  +  s^^u  )2it  r  Az  At  . 

\  r  z  rr  /  1 


This  completes  the  discussion  of  the  Phase  3  calculation 
for  a  typical  cell,  of  the  momentum  and  total  energy  increments 
due  to  the  deviator  stresses  acting  on  the  cell.  The  two  momen 
turn  increments  are  used  to  determine  the  velocity  increments 
by  dividing  by  cell  mass.  This  fixes  the  Phase  3  change  in 
cell  kinetic  energy  j  that  the  internal  energy 

can  be  calculated  using  the  updated  value  of  total  energy 
Err  =  E  +  i  u-u..  At  the  end  of  Phase  3  there  are  new  values 
of  E^Ck),  u(k),  v(k)  as  updated  by  Phase  3.  As  in  Phase  1, 
the  change  in  internal  energy  for  mixed  cells  is  partitioned 
among  the  materials  in  proportion  to  their  partial  volumes. 


2.2.4  Phase  2 — The  Effects  of  Transport 

The  purpose  of  this  section  is  to  describe  how  the 
transport  of  mass,  momentum  and  energy  from  cell  to  cell  is 
accounted  for  in  the  code.  This  is  done  by  calculating  the 
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integrals  in  the  last  terms  o£  Eqs.  (10)  through  (12).  In 
the  discussion  below,  primary  attention  is  given  first  to 
the  case  of  transport  between  pure  cells.  At  the  conclusion 
of  this  discussion,  the  special  provisions  required  to  treat 
mixed  cells  are  described. 

2 . 2 . 4 . 1  Continuity  Equation,  (10)  ^ 

The  transport  mass  is 

I 

(m)  =  -  At  J" pu^n^ds 

s  i 

and  is  determined  for  each  cell  face  from  ! 

6m  =  -  u^A^At  . 

is  the  density  of  the  cell  from  which  the  mass  moves 
(donor  cell) ,  A^  is  the  area  of  the  face  and  u^  is  an 
interpolated  value  of  the  velocity  component  normal  to  A^ 
representing  approximately  the  velocity  at  the  interface 
at  the  end  of  the  time  step.  For  example,  considering 
cells  k  and  ka  above, 

^[v"ck)  t  v"Cka)]  ; 

~  [i  .  vl(k)  ’ 

Calculated  transport  masses  are  subtracted  from  the  donor 
cell  mass  and  added  to  the  acceptor  cell  mass.  This  updating 
is  done  after  the  transport  terms  have  been  calculated;  all 
transport  terms  are  computed  using  the  pre- transport  quantities. 
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2. 2.4. 2  Equation  of  Motion,  CH) 

1.  Axial  Motion:  The  term  in  Eq.  (11)  for  axial 
momentum  transport  is 

A^(mu.)  =  -  At  y*  (pu.u.)n.ds  . 
s 

At  each  face  of  the  cell  the  transport  specific  momen¬ 
tum,  Uj ,  is  taken  to  be  the  axial  velocity  of  the  cell  from 
which  the  mass  moves  (donor  cell,  index  kd) ,  i.e., 

u^  =  v^(kd)  . 

Since  the  various  faces  of  the  cell  have  different  donor  cells 
it  is  convenient  to  express  the  momentum  transport  for  each 
face 


6 (mv)  =  v^(kd)6m 


where 


6m  =  -  u . A . At 

is  the  mass  which  is  transported  across  the  cell  face,  as 
given  in  the  preceding  mass  transport  calculation.  Note  that 
5(mv)  is  the  axial  momentum  and  that  6m  is  the  mass  transport 
in  either  the  r  or  z  direction,  depending  on  which  face  of 
the  cell  is  being  computed. 

2.  Radial  Motion:  Again  from  Eq.  (11), 

A^(mu.)  =  -  At/* (pu.u.)n.ds 

•^s 
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and,  by  analogy  with  the  axial  case,  the  equation  for  the 
transport  of  radial  motion  across  an  interface  is 

6  (mu)  =  u’^(kd)6m 

where  u^(kd)  is  the  time  n  velocity  of  the  donor  cell  and 
6m  =  -p^  u^A^At  is  the  mass  which  is  transported  across  the 
cell  face  in  question,  as  computed  in  the  continuity  equation 
(above) . 

2 . 2 . 4 . 3  Energy  Equation  (12) 

The  Eq.  (12)  expression  for  the  transport  of  energy 
is 


To  evaluate  this  integral,  the  transported  specific  energy  is 
taken  to  be  that  of  the  donor  cell,  kd,  i.e.. 


and  the  total  energy  which  is  transported  across  a  given 
interface  is  therefore  the  product  of  this  specific  energy 
and  the  associated  transport  mass  which  was  computed  above 
in  the  continuity  equation,  i.e., 

AalmEf)  =  Ep6m  . 

Once  the  mass,  momentum  and  energy  transports  are  known 
for  all  faces  of  the  cell,  the  cell  quantities  can  be  updated 
for  these  effects.  New  cell  velocities  are  determined  by 
dividing  the  new  momenta  by  updated  cell  mass.  This  fixes 


i 

] 
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the  new  kinetic  energy  so  that  the  internal  energy  can  be 
calculated  from  the  known  total  energy. 

2. 2. 4. 4  Extension  of  Phase  2  to  Mixed  Cells 

The  procedure  used  to  treat  mixed  cells  is  an  extension 
of  the  preceding  method  for  transport  between  pure  cells. 
Specifically,  the  mass,  say  for  material  b,  which  is  trans¬ 
ported  across  a  given  interface  is,  in  analogy  with  the  re¬ 
sult  in  Section  2. 2. 4.1  above^  just 

s 

where  now  is  the  donor  cell  density  for  material  b  and 

the  surface  integral  is  now  over  just  that  part  of  the  surface 
area  belonging  to  material  b.  The  determination  of  the 
density  is.  made  in  the  pressure  iteration,  to  be  the 

subject  of  Section  V;  the  determination  of  the  fractional 
surface  area  belonging  to  a  given  material  is  part  of  the 
discussion  in  Section  IV  of  the  method  for  computing  the 
intersection  of  material  interfaces  and  cell  boundaries. 

The  momentum  transport,  in  analogy  with  the  pure  cell 
case,  is  just  this  mass  transport  multiplied  by  the  appro¬ 
priate  velocity  component  (axial  or  radial)  of  the  donor  cell 
material. 

The  energy  transport  is  the  mass  transport  multiplied 
by  the  energy  per  unit  mass,  kinetic  plus  internal.  The 
kinetic  energy  per  gram  is  the  same  for  all  materials  in  the 
donor  cell,  but  the  internal  energies  are  generally  different 
(internal  energies  were  updated  in  Phase  1  and  Phase  3  as 
described  in  prior  sections) . 
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As  in  the  case  of  pure  cells,  final  values  of  cell 
variables  are  easily  determined  once  all  the  transport  quan¬ 
tities  are  known  for  the  cell.  The  new  cell ^velocity  com¬ 
ponents  are  determined  by  dividing  the  new  momenta  by  the 
updated  cell  mass.  This  fixes  the  new  kinetic  energy  for 
each  material;  the  internal  energy  for  each  material  is  then 
chosen  to  conserve  total  energy,  kinetic  plus  internal,  for 
the  individual  materials. 


2.3  THE  EQUATION  OF  STATE 
2.3.1  The  Determination  of  P 


The  equation  of  state  used  in  HELP,  for  initially  con¬ 
densed  materials,  is  that  due  to  J.  H.  Tillotson, modified 
to  give  a  smooth  transition  between  condensed  and  expanded 

states.  For  the  condensed  states,  i.e.,  when  p/p  >  1, 

0 

or  for  any  cold  states,  I  the  equation  has  the  form 


a '  + 


Ip 


Au  +■ 


For  expanded  hot  states,  i.e.,  when  p/p  <  1:  and  I  > 

0  s 

the  equation  of  state  has  the  form 


P  =  Pg  =  alp  + 


bip 


LI 

0. 


+  1 


+  Aye 


-6(V/V  -l) 


■a(v/V  -if 


A  smooth  transition  between  the  condensed  and  expanded  states 
is  insured  by  a  transition  equation  for  the  intermediate  region 
defined  by  I  <  I  <  and  p/p  <1.  This  blended  portion 

S  S  Q 
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o£  the  equation  o£  state  has  the  £orm 

K-h 

In  these  equations,  p,  I  and  p  are  pressure,  speci£ic  inter¬ 
nal  energy  and  mass  density  respectively,  ri  =  p/p  =y  +  l  =  V  /V, 

0  0 

and  p  ,  a,  b,  I  ,  I  ,  1'  A,  B,  p  and  p  are  constants  £or 
0  0  ^  is  0 

the  particular  material.  The  values  of  these  constants  for 
nineteen  materials  are  given  in  Table  I. 

For  mixed  cells,  additional  considerations  are'  require 
to  determine  the  cell  pressure  and  the  densities  of  the 
various  materials  within  the  cell.  This  problem  is  discussed 
as  Section  V. 

2.3.2  The  Determination  of  Deviatoric  Stresses  . 

The  deviatoric  stress  increments,  ds^j ,  are  determined 
by  using  the  elastic  relation 

ds.  .  =  2Gder  . 

where  G  is  the  modulus  of  rigidity  and  are  the  incre¬ 

ments  of  deviatoric  strain.  When  such  an  increment  of  stress 
causes  the  yield  criterion  to  be  violated  each  stress  component 
is  proportionately  reduced  to  bring  the  stress  state  normally 
back  to  the  yield  surface.  A  variable  yield  strength 

Y  =  (k^>  K^p  .  K^p")(l  -  E/E„) 

is  defined,  to  account  for  the  increase  in  strength  at  high 
pressures  and  the  decrease  in  strength  at  elevated  temperatures. 
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TABLE  I 

HELP  EQUATION  OF  STATE  CONSTANTS^ 


Material 

Code 

Number 

Material 

a 

b 

A 

2 

dynes/cm 

B 

dynes/cm^ 

^0 

ergs/g 

a 

8 

Es 

P 

g/cm^ 

— 

1 

w 

.5 

1.04 

3.08  w  10^^ 

2.5  X  10^^ 

.225  X  10^^ 

10 

10 

1.11  X  10^° 

19.17 

5.6  X  10^® 

2 

Cu 

.5 

1.5 

1.39 

1.1 

.325 

5 

5 

1.38 

8.9 

6.9 

3 

Fe 

.5 

1.5 

1.28 

1.05 

.095 

5 

5 

2.44 

7.9 

10.2 

4 

Al 

.5 

1.63 

.75 

.65 

.05 

5 

5 

3.0 

2.7 

15.0 

5 

Be 

.55 

.62 

1.17 

,55 

.175 

5 

5 

10,0 

1.8 

46.0 

6 

Ti 

.5 

.60 

1.03 

.5 

.07 

5 

5 

3.5 

4.5 

12.5 

7 

Ni 

.5 

1.33 

1.91 

1.5 

.09 

5 

5 

2.85 

8.9 

9.4 

8 

Mo 

.5 

1.02 

2.71 

1.65 

.045 

5 

5 

2.8 

10. '2 

9.0 

9 

Th 

.4 

.86 

.53 

.5 

.025 

9 

.88 

2.0 

11.7 

8,0 

10 

Pb 

.4 

2.4 

.466 

.15 

.02 

10 

2 

.26 

11.3 

.97  X  10'* 

11 

CH^ 

.6 

2.0 

.075 

.02 

.07 

10 

5 

2.4 

.9 

18.0 

12 

2 

Granite 

.5 

1.3 

.60 

.00 

.16 

5 

5 

3.5 

2.7 

.18 

13 

Andesite^ 

.5 

1.3 

.34 

.28 

.16 

5 

5 

3.5 

2.7 

.18 

14 

Wet  Tuff 

.5 

1.3 

.10 

.06 

.11 

5 

5 

3.2 

.  2.0 

.16 

15 

Dry  Tuff 

.5 

1.3 

.045 

.03 

.06 

5 

5 

3.2 

1.7 

.18 

16 

Oil  Shale 

.  5 

1.0 

2.8,  _ 

.11 

.11 

_,5 

-A  - 

3.2 . . 

...  2A  ,  , 

.16 

17 

Dolomite 

.5 

.6 

.85 

,30 

.10 

5 

5 

2.5 

2.8 

.14 

18 

Limestone 

.5 

.6 

,4. 

.67 

.10 

5 

5 

2.5 

2.7 

19 

Malite 

.5 

.6 

.25 

.30 

.05' 

5 

5 

2.0 

2.2 

.15 

^Dienes,  J.  K.,  et  al . ,  "An  Eulerian  Method  for  Calculating  Strength  Dependent  Deformation,*' 

General  Atomic  Report  GAMD-8497,  February  1968, 

2 

For  these  two  materials  the  constants  are  valid  only  for  pressures  less  than  150  kbar. 

For  a  formulation  of  the  high  pressure  equation  of  state  of  these  materials,  see 

Allen,  R.  T.,  "Equation  of  State  of  Rocks  and  Minerals,”  General  Atomic  Report  GAMD-7834,  March  1967. 
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2.3.3  Material  Failure  in  Tension 

The  failure  criterion  is  a  simple  one  based  upon 
relative  volume.  When  the  relative  volume  in  a  cell  reaches 
a  certain  value  which  is  greater  than  a  specified  maximum 
distension,  that  cell  is  said  to  fail  and  any  computed  ten¬ 
sions  are  zeroed  out. 

2.4  OTHER  FEATURES  OF  THE  NUMERICAL  METHOD 

2.4.1  Grid  Boundaries  and  Axis 

Additional  specification  is  required  for  cells  which 
border  the  grid  boundaries  or  axis  of  symmetry,  because 
necessary  quantities  are  not  defined  for  neighbor  cells  which 
would  be  outside  the  grid.  The  code  provides  for  a  trans- 
mittive  boundary  at  the  top  and  right,  an  option  of  a  trans- 
mittive  or  reflective  boundary  at  the  bottom  and  a  reflective 
boundary  at  the  axis.  Boundary  conditions  for  border  cells 
are  then  derived  by  assuming  (fictional)  neighbor  cells  out¬ 
side  the  grid.  For  transmittive  boundaries  the  flow  variable 
are  the  same  in  the  fictional  cell  as  in  the  border  cell,  and 
for  reflective  boundaries  the  state  is  assumed  to  be  the  same 
except  that  the  velocity  component  normal  to  the  boundary  has 
the  opposite  sign. 

2.4.2  Re  zone 

A  rezone  subroutine  can  be  triggered  by  mass  motion  out 
of  the  right  or  top  grid  boundaries.  Four  cells  are  combined 
into  one,  starting  at  the  lower  left  corner  (z  =  0,  r  =  0)  of 
the  grid  and  new  cells  are  added  at  the  top  and  right.  The 
total  number  of  cells  and  their  shape  is  held  constant.  In 
combining  cells  during  rezone,  total  mass,  momentum  and  energy 
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are  conserved.  The  net  effect  is  a  slight  increase  in  inter¬ 
nal  energy,  corresponding  to  the  kinetic  energy  which  the  four 
cells  had  in  a  coordinate  system  moving  with^  their  center  of 
mass. 

In  some  cases,  and  specifically  when  a  shock  passes 
through  a  transmittive  grid  boundary,  it  has-been  found  ; 
necessary  to  rezone  (thus  keeping  the  entire^  disturbance 
within  the  new  grid)  to  avoid  undesirable  oscillations  in 
the  calculation.  The  oscillations  apparently  indicate  that 
a  simple  transmittive  boundary  condition  does  not  adequately 
simulate  the  effect  of  additional  material. 

2.4.3  Variable  Zoning 

I 

t 

Although  the  discussion  proceeded  with  no  special 
reference  to  variable  zoning,  the  code  has  the  capability 
of  computing  problems  in  which  the  zone  size  is  variable. 

If  constant  zones  are  not  specified,  the  cell  sizes  must 

be  supplied  as  a  part  of  the  SETUP  deck.  It  is  possible  to 

call  for  variable  cell  size  in  one  direction  : and  constant 
in  the  other. 

The  zones  should  be  selected  so  that  the  discontinuities 
in  size  between  adjacent  cells  are  not  excessive  and  the  aspect 
ratio  of  the  cells  are  not  too  extreme.  It  was  found  that  an 
initially  spherical  explosion  loses  its  symmetry  in  the  course 
of  a  calculation  when  the  ratio  of  the  cell  edges  is  greater 
than  four  to  one,  and  that  satisfactory  symmetry  is  obtained 

when  the  ratio  is  two  to  one.  Cell  size  variations  are 

usually  held  to  less  than  10  percent  between  adjacent  cells. 
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III.  NUMERICAL  METHOD  FOR  ELASTIC-PLASTIC  EFFECTS 


3.1  INTRODUCTORY  REMARKS 

The  effects  of  material  strength  are  computed  as 
Phase  3,  and  when  this  phase  is  simply  by-passed  the  remaining 
code  solves  problems  in  compressible  fluid  dynamics.  Phase  3 
adds  the  resistance  to  shear  which  characterizes  solids. 

The  Phase  3  determination  of  strength  effects  is  it¬ 
self  divided  into  two  principal  parts.  One  of  these,  the 
effects  of  (predetermined)  deviatoric  stresses  in  changing 
the  material  velocity  field  and  the  material  specific  inter¬ 
nal  energies,  has  been  described  above  as  Section  2.2.3, 

The  second  part,  discussed  below,  is  the  determination  of 
these  stress  deviators,  given  the  velocity  field  and  the 
material  constitutive  equations. 


3.2  DETERMINATION  OF  STRAIN  RATE  DEVIATORS 

As  a  first  step  in  computing  stress  deviators,  the 
instantaneous  strain-rate  deviators  are  computed  from  the 
velocity  field,  i.e.,  for  cylindrical  coordinates 


e  =  ( u  +  V  )/2 

rz  V  y  xr 

Here  the  necessary  space  derivatives  of  the  velocity  field 
are  computed  (centered  at  cell  centers)  in  a  straightforward 
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I  ■  “  “1 

manner ,  i . e . ,  I  I 


I  kb  1 
1  I 

I _ __1 


But  the  above  calculation  gives  only  the  strain  rate  devia¬ 
tor  for  that  material  which  is  at  cell  center  at  the  beginning 
of  the  time  step,  and  it  is  necessary  to  recognize  that  con¬ 
stitutive  equations  are  to  be  applied  only  along  a  particle 
path.  To  this  end,  one  must  take  a  convective  derivative. 
Consider  the  time  N  position  of  a  particle  which,  at 


time  N+1,  will  be  at  cell  center.  This  point  will  be  offset 


Now  we  wish  to  determine.,  by  interpolation  using  the  value  of 
the  strain  rate  deviators  at  cell  centers,  the  strain  rate 
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deviators  at  the  offset  point  P^.  This  is  done  by  weighting 

the  contribution  of  neighboring  cells  in  proportion  to  their 

overlap  areas  with  a  rectangle  of  cell  dimensions,  centered 
N 

at  P  ,  For  example,  in  the  above  sketch,  the  weighting 
factors  are 

for  cell  k 

WK  =  (Ax  -  1 6x1) (Ay  -  l6y|) 
for  cell  in  horizontal  direction 
WKH  =  I 6x| (Ay  -  I 6y I ) 
for  cell  in  vertical  direction 
WKV  =  I 6y| (Ax  -  I 6x| ) 

for  cell  in  diagonal  direction  •  ' 

WKD  =  I  6x1  I Syl 

so  that  the  resulting  interpolated  value  of  a  strain  rate 
deviator  component,  say  given  by 

(WK)£^^(k)  +  (WKH)c^Jkh)  +  (WKV)  (kv)  +  (WKD)  (kd) 
Sz  WK  +  WKH  +  WKV  +  WKD 

Having  thus  determined  e  ,  e  ,  e  for  the  mass  which  will 
^  zz’  rr’  rz 

be  at  the  center  of  cell  k  at  the  end  of  the  time  step,  we 
can  now  update  the  strain  deviators.  Say,  for 
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3.3  DETERMINATION  OF  STRESS  DEVIATORS 


These  strain  increments  are  next  used  in  the  material 
constitutive  equation,  to  update  the  stress  deviator  tensor 
s^j.  To  do  this  one  needs  the  time  N  stress  deviator  for 
the  particle  in  question  and  these  are  determined  (from  stored 
cell  centered  values  of  the  at  time  N)  by  an  identical 

area  weighting  procedure  as  described  above  for  strain  rate 
deviators . 

For  an  elastic-plastic  material  the  following  calculation 


then  (provisionally)  updates  the 


according  to 


N+1 


=  s 


N 

ij 


+  2G  Ae.  - 
iJ 


where  G  is  the  material  rigidity  modulus,  i.e.,  the  three 
components  of  interest  are  given  by 


N+1 

s 

Z2 


s^  +  2G  Ae 
z  z  zz 


^N+1 

s 

rr 


2G  Ae 


rr 


sN+1 

^r  z 


2G  Ae 


r  z 


The  provision  (above)  is  that  the  yield  condition  for  the 
material  not  be  violated,  i.e.,  one  tests  to  see  if  the 
second  invariant  exceeds  the  yield  condition 


s  .  .  s  .  . 
11  iJ 


=  S' 


zz 


+  s' 


rr 


+  s ' 


66 


2s^ 
r  z 


>  i2Y‘ 


Here  S-_  =  - (s  +  s  )  and  Y  is  the  yield  strength  in 
yy  VTTzz/ 

shear.  As  pointed  out  in  Section  2.3,  the  yield  point  can 
be  some  function  of  the  local  thermodynamic  state  (i.e.,  of 


.1 
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p  and  E) .  I£  the  yield  point  is  exceeded  during  a  time 
step,  as  would  occur  in  plastic  deformation,  all  of  the 
stresses  are  proportionately  reduced  to  cause  the  yield  con¬ 
dition  to  be  satisfied.  The  second  invariant  being  smaller 
than  2Y^,  on  the  other  hand,  corresponds  to  elastic  deforma¬ 
tion,  for  which  the  stresses  are  not  reduced. 

The  preceding  calculation  updates  the  s^^ ,  giving 
values  of  these  quantities  at  cell  centers  (as  corrected  for 
material  convection.)  For  the  use  of  stress  deviators  in 
updating  velocities  and  energies,  however,  the  method  described 
as  Section  2.2.3  assumes  values  of  these  quantities  at  cell 
surfaces  rather  than  cell  centers.  These  cell-surface  quan¬ 
tities  are  obtained  by  merely  averaging  the  corresponding 
quantities  at  the  centers  of  the  two  cells  in  question. 

3.4  MODIFICATIONS  FOR  GRID  BOUNDARIES  AND  AXIS 

The  description-  given  above  applies  to  normal  cells 
within  the  grid.  Additional  statements  must  be  made  to  cover 
grid  boundaries  and  the  axis  of  symmetry. 

Grid  boundaries  may  be  reflective  (bottom  grid  boundary) 
or  transmittive  (bottom,  right  or  top) .  In  calculating  the 
strain  rate  deviators  for  reflective  cells  (from  the  velocity 
field)  it  is  assumed  that  the  cell  which  would  be  outside  the 
grid,  has  the  same  tangential  component  of  velocity  and  an 
equal  and  opposite  normal  component  of  velocity,  as  the  cell 
within  the  grid.  Thus,' for  cell  k  at  a  reflective  bottom 
grid  boundary,  the  assumed  cell  just  below  k  has  velocities 

u  =  u(k) 

V  =  -v(k) 

These  velocities  are  then  used  in  an  otherwise  standard  (above) 
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calculation  of  the  strain  rate  deviator  tensor  for  cell  k.; 

For  transmittive  boundary  cells,  one  assumes, that  the  cell 
outside  the  grid  has  both  of  its  velocity  components  equal 
to  the  velocity  of  the  cell  just  inside  the  boundary. 

u  =  u(k) 

v  =  .v(k) 

With  these  assumptions,  it  is  then  possible  to  compute  strain 

i 

rate  deviators  for  the  border  cells.  For  either  transmittive 
or  reflective  grid  boundaries,  the  normal  stress  at  the  grid 
boundary  is  taken  to  be  that  of  the  border  cell.  The  tangen¬ 
tial  stress  is  zero  for  reflective  boundaries  and  is  equal 
to  that  in  the  border  cell  for  transmittive  boundaries;  e.g., 
for  a  reflective  boundary  at  grid  bottom  ■ 

J 

SNB  =  s^^(k)  stress  normal  bottom 

STB  =  0  stress  tangent  bottom 

and  for  a  transmittive  boundary  at  the  right  of  the  grid 
SNR  =  s^^(k)  stress  normal  right 

STR  =  s^^(k)  stress  tangent  right 

The  axis  of  symmetry  is  essentially  a  reflective  grid 
boundary.  In  calculating  strain  rate  deviators  for  the  column 
of  cells  next  to  the  axis,  the  cell  "to  the  left"  is  assumed 
to  have  the  same  axial  velocity  and  an  equal  and  opposite 
radial  velocity,  i.e., 

V  =  vCk) 

u  =  -u(k)  ' 
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Stresses  for  the  axial  boundary  of  this  column  of  cells  are 
simply  set  to  zero  in  the  program,  since  the  cell  area  over 
which  these  stresses  act  is  zero  and  it  is  the  product  of 
stresses  and  areas  which  is  of  importance  in  computing  the 
effects  on  cell  velocities  and  energy  (Section  2.2.3). 

3.5  MODIFICATIONS  FOR  MATERIAL  FREE  SURFACES 


Cells  within  the  grid  which  contain  a  material  free 
surface  require  special  consideration.  In  general,  such 
cells  are  only  partly  filled  with  material  and  are  partly 
void.  The  program  complications  that  would  be  involved  by 
detailed  treatment  of  such  partially  filled  cells  are  avoided 
by  simply  insuring  that  such  cells  have  the  same  velocity 
and  specific  internal  energy  as  their  highest  pressure 
neighbor  cell.  This  is  done  in  the  GLUE  routine  (see 
Section  7.1.10),  after  Phase  1  is  executed.  The  glueing  of 
two  such  cells  is  done  by  conserving  total  momentum  and  energy 
of  the  pair  of  cells,  and  leaving  the  mass  within  each  cell 
the  same.  After  glueing,  the  two  cells  have  the  same  velo¬ 
city  and  specific  internal  energy.  In  Phase  3,  then,  the 
stress  deviators  are  assumed  to  be  zero  within  a  partially 
filled  free  surface  cell  and  the  stresses  at  the  boundaries 
between  this  cell  and  its  neighbors  are  also  assumed  to  be 
zero.  The  combined  effect  of  glueing  and  neglecting  actual 
stresses  would  be  sufficient  to  cause  the  free  surface  cell 
material  to  have  the  same  velocity  as  the  main  body  of 
material . 
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IV.  THE  PROPAGATION  OF  TRACER  PARTICLES 
AND  THEIR  USE  TO  DESCRIBE  MATERIAL  INTERFACES 

4.1  INTRODUCTORY  REMARKS 

The  most  complicated  logic  of  the  HELP  code  is  found 
in  the  subroutines^  dealing  with  the  definition  and  use  of 
material  interfaces.  Conceptually  the  task  is  not  difficult. 

The  idea  of  moving  tracer  particles  with  a  locally  defined 
velocity  and  of  using  them  to  describe  the  interface  between 
materials  is  a  simple  one.  Likewise,  employing  the  partial 
cell  boundary  areas  determined  by  interface  position  to  com¬ 
pute  the  mass  flux  of  materials  is  a  very  natural  use  of  these 
interfaces.  The  complexity  of  the  task  comes  with  making  the 
treatment  completely  general.  As  the  code  is  written,  there 
are  no  restrictions  on  the  direction  of  the  flow,  on  the 
shape  of  the  material  interface,  on  the  number  of  times  an 
interface  can  cross  a  cell  boundary,  or  on  the  number  of 
materials  contained  in  a  cell. 

It  is  hoped  that  the  present  section  will  aid  the  user 
in  understanding  the  organization  and  logic  of  these  routines. 

In  particular,  this  section  covers  the  definition  of  the 
material  interfaces,  the  computation  of  fractional  cell  boun¬ 
dary  areas  for  each  material  in  a  mixed  cell,  and  the  use  of 
these  fractional  areas  in  defining  the  mass  transport  terms. 

The  creation  of  mixed  and  pure  cells,  and  the  propagation  of 

2 

the  tracer  particles  that  define  the  package  boundaries  will 
also  be  discussed. 

^INFACE,  NEWMIX,  NEWRHO,  FLUX. 

^To  distinguish  the  INFACE  and  PH2  [Phase  2)  functions,  it  should 
be  noted  that  the  mass  transports  (for  mixed  cells)  are  only 
computed  and  stored  in  INFACE.  The  actual  transport  of  all 
quantities,  and  the  associated  updating  of  cell  variables,  is 
performed  in  PH2  for  both  mixed  and  pure  cells. 
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4.2  DEFINING  MATERIAL  INTERFACES 

A  set  of  massless  tracer  particles  are  initially  posi¬ 
tioned  along  the  boundary  of  each  material  package.  These 
tracers  are  numbered  such  that  the  package  is  on  the  left  as 
one  proceeds  between  any  two  consecutive  tracers.  The  material 
interface  therefore  is  defined  by  the  line  segments  connecting 
these  tracers.  The  initial  density  of  tracers  is  controlled 
by  input  variables  and  can  be  increased  in  a  specified  region 
as  the  calculation  progresses.  (See  description  of  ADDTCR, 
Section  7.1.1.) 

To  prevent  the  boundaries  of  adjacent  packages  from 
crossing,  the  tracers  along  the  common  boundary  coincide 
exactly.  Because  these  identical  points  are  moved  with 
exactly  the  same  velocity,  they  remain  superimposed  diiring 
the  entire  calculation. 

To  allow  for  spatially  disconnected  subpackages,  the 
code  employs  the  following  convention.  The  last  point  of  a 
package  or  subpackage  is  a  dummy  point  and  is  identical  to 
the  first  point  of  the  package  or  subpackage.  The  line 
drawn  between  the  last  two  points  is  therefore  not  considered 
part  of  the  package  boundary.  Likewise  that  ;last  point  is 
assumed  to  be  disconnected  from  the  first  point  of  the  next 
subpackage  if  such  a  subpackage  exists. 
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4.3  COMPUTING  FRACTIONAL  AREAS 
4.3.1  Problem  Definition 


When  an  Eulerian  cell  is  cut  by  the  material  N 
interface,  we  are  interested  in  knowing  what  fraction  of 
each  side  of  the  cell  should  be  associated  with  material 

N. 


For  example,  the  fractional  areas  associated  with  material 
N  in  the  diagram  above  are  as  follows: 

Left  =0. 

Top  .8  TT^X?  - 

Right  .7(y^  -  yj.i)  2tt  x^ 

Bottom  .3ir  |  x?  - 

These  fractional  areas  are  used  to  compute  the  mass  flux  of 
material  N  across  each  side  of  the  mixed  cell.  (Actually, 
for  each  mixed  cell  it  is  necessary  to  compute  and  store  only 
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■| 

t 

I 

the  fractional  areas  for  the  right  and  top  sides  since  the 
cells  below  and  on  the  left  provide  the  information  for  the 
other  two  sides.)  The  method  given  below  for  computing 
these  fraction-areas  puts  no  restrictions  on!  how  many  inter¬ 
faces  are  in  a  cell  or  on  how  many  times  a  single  interface 

1 

crosses  a  cell  or  cell  boundary.  It  does,  however,  assume 
that  a  material  package  boundary  never  crosses  itself,  and 
experience  has  shown  that  when  it  does  the  logic  of  INFACE 
breaks  down.  This  is  not  a  limitation  of  the  method  since 
the  boundary  will  not  cross  itself  unless  the  tracers  become 
very  sparse,  and  this  can  be  prevented  by  adding  new  tracers, 
as  discussed  in  the  description  of  the  subroutine  ADDTCR. 

4.3.2  Computing  Intersections  of  Interfaces! 

witn  Grid  Lines  : 

■  ■  ■  “  I 

I 

Subroutine  INFACE  considers  consecutively  (two  at  a 
time)  the  tracer  points  which  circumscribe  each  material 
package  in  a  counter-clockwise  direction.  For  simplicity, 
the  coordinates  of  the  tracers  are  in  grid  line  units  rather 
than  centimeter  units.  Thus  the  coordinates ! of  a  pair  of 
tracers  immediately  indicate  whether  they  straddle  one  or 
more  cell  boundaries. 

TX  and  TY  are  the  FORTRAN  variable  names  of  a  tracer’s 
X  and  y  coordinates,  respectively.  The  arrays  are  doubly 
dimensioned;  the  first  dimension  identifies  the  material 
package,  the  second  gives  the  sequential  ordering  of  the 
points.  In  the  example  shown  below,  the  indices  indicate  that 
both  tracer  points  belong  to  material  package  one  (first 
index)  and  that  they  are  consecutive  points  (second  index) , 
namely  the  20th  and  21st  points  describing  package  one.  The  x- 
coordinates  indicate  the  line  through  the  points  does  not  inter¬ 
sect  a  vertical  grid  line,  whereas  the  y-coordinates  indicate  the 
line  through  the  points  intersects  the  fourth  horizontal  grid 
line. 

i 
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TXC1,20)  =  2.5 
TY(1,20)  =  3.5 
TX(1,21)  =  2.8 
TY(1,21)  =  4.6 


4.3.3  Defining  Fractional  Areas  of 

Intersected  Cell  Boundaries 

To  continue  with  the  above  example,  INFACE  defines  the 
point  of  intersection  by  a  straight  line  interpolation  between 
tracers,  and  uses  it  to  define  the  fractional  area  of  the . top 
cell  boundary  associated  with  material  package  one.  The 
code  recognizes  from  the  sequential  numbering  of  the  points 
that  material  one  is  on  the  left  of  the  point  of  intersection. 

The  FORTRAN  variables  FRACTP  and  FRACRT  store  the 
fraction  areas  of  the  top  and  right  boundaries,  respectively 
for  all  the  materials  of  a  mixed  cell.  These  arrays  are  also 
doubly  dimensioned.  Given  FRACTP (N ,M)  ,.  FRACRT(N,M),  the  first 
dimension,  N,  identifies  the  material  package  number,  and  the 
second  dimension,  M,  is  an  index  that  links  the  mixed  cell  K 
to  the  special  mixed  cell  arrays  1^7  the  relation,  M=FLAG  (K) -100 . 
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If  material  package  two  is  adjacent  to  material  package 
one,  the  material  two  tracers  are  identical  to  the  material 
one  tracers  but  are  in  an  opposite  order.  For  example,  package 
two  tracers  (say  30  and  31)  that  coincide  with  the  package  one 
tracers  discussed  previously  would  be  as  follows: 

TX(2,30)  =2.8 

TY(2,30)  =4.6 


TYC2,31)  =  2.5 
TYC2,31)  =3.5  1 

The  fractional  area  to  the  right  of  the  intersection  point  is 
associated  with  material  package  two,  and  is  stored  in 
FRACTP(2,M) . 

4.3.4  Defining  Fractional  Areas  of  Boundaries 
Not  Intersected" 

Because  the  fractional  area  is  a  term  in  the  flux  equa¬ 
tions  applied  to  mixed  cells,  it  must  be  defined  for  the  right 
and  top  boundary  of  every  mixed  cell  even  if  one  or  both  are 
not  intersected  by  a  material  interface.  In  that  case  the 
fractional  area  either  equals  the  total  cell  boundary  area 
or  is  zero  ,  e . g .  ! 


?1aterial  Package  1  FRACTP(1,M)  =  jr  *  (x?  -  x?_^J 


Material 
Package  2 


\ 


FRACRT(1,M)  =  (y^ 

FRACTPC2,M)  =  :0. 
FRACRTC2,M)  =  0. 


X  . 

1 
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In  the  above  example,  none  o£  the  materia]  of  package  two  is 
transported  across  the  right  or  top  boundaries  of  the  cell. 

On  the  other  hand,  the  full  cell  area  is  used  to  compute 
the  mass  flux  of  the  material  of  package  one  across  these 
boundaries.  There  is  clearly  a  need  to  systematically  define 
the  fractional  areas  of  mixed  cell  boundaries  that  are  not 
intersected  by  a  material  interface.  The  following  discussion 
illustrates  how  the  fractional  area  of  these  non-intersected 
boundaries  are  correctly  defined  by  "presetting"  them  when  the 
interface  enters  a  cell  and  "resetting"  them  when  it  leaves. 

The  direction  of  an  interface,  given  by  the  sequential 
order  of  the  tracer  points,  determines  whether  one  is  entering 
a  cell  or  leaving  it.  When  entering  a  cell,  and  if  crossing 
this  cell  boundary  for  the  first  time,  the  program  processes 
tl\e  other  sides  of  the  cell  in  a  clockwise  order  and  presents 
their  fractional  areas  to  equal  the  total  cell  boundary  area 
until  it  comes  to  a  side  that  has  already  been  crossed  by  the 
same  material  interface  (i.e.,  a  side  which  has  an  associated 
fractional  area  that  is  non-zero  and  is  less  than  the  total 
area. ) 

Fractional  areas  preset: 

Bottom  =  "  ^i  l) 

Left  =  (y.  -  yj.i)2Tr  x-.^ 

Top  =  ^(x?  -  x?.J 
i 


Case  (1) 


i-1 
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Fractional  areas  preset: 
Bottom  = 


I£  the  interface  of  material  N  crosses  a  cell 
boundary  for  the  second  time,  the  code  does  not  preset 
of  fractional  areas  when  material  N  is  between  the  two 
points  of  intersection  (Case  4) .  If  material  -  N  is  out¬ 
side  the  two  points,  the  fractional  areas  are  preset  as 
described  above  (Case  5),  The  program  senses  that  material 
N  is  between  the  points  when  the  sum  of  the  fractional  area 
computed  on  the  previous  crossing (s)  and  the  one  computed 
currently  is  greater  than  the  total  area. 
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Actually  this  rule  holds  regardless  of  how  many  times  the 
interface  crosses  the  boundary,  as  illustrated ‘by  Case  6  and 
Case  7.  i 


Case  (6)  \ 

j  r- 


Previous 


revious 


Current 


of  areas  <  Total  area 


Fractional  areas  preset: 
Bottom  = 

Left  =  (y.  -  X..J 


Case  (7) 


1  ~~k 


^  ^^revi 


Current 


Sum  of  areas  >  Total  area 


Fractional  areas  preset: 


None 


When  the  interface  enters  cell  K  from  the  right, 
it  has  just  left  cell  K+1,  and  the  program  must  "reset"  to 
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zero  the  fractional  areas  of  certain  non-intersected  sides 
of  cell  K+1. 


In  cases  1  through  3  the  side  where  the  interface  leaves 
is  crossed  only  once  by  this  interface.  Proceeding  in  a  clock¬ 
wise  direction  from  the  side  the  interface  has  left,  the  program 
resets  to  zero  the  fractional  areas  of  the  non-intersected  sides 
until  it  encounters  a  side  that  is  intersected  by  this  interface. 
Case  (1) 


Fractional  areas  reset: 

Top  =0. 

Right  =  0. 


Case  (3) 


Fractional  areas  reset: 
None 


In  cases  4  and  5  below,  the  side  where  the  interface 
leaves  has  been  previously  crossed  by  this  interface.  The 
program  sums  the  fractional  areas  resulting  from  a  previous 
crossing  of  that  boundary  and  the  current  crossing.  If  the 
sum  of  these  areas  is  less  than  the  total  cell  boundary  area 
(Case  4)  none  of  the  fractional  areas  of  the  other  sides  are 
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reset.  I£  the  sum  is  greater  (Case  5),  the  program  proceeds  to 
reset  to  zero  the  areas  o£  the  other  non-intersected  sides  o£ 


<  Total  area 


None 


Case  (5) 


areas  >  Total  area 


Fractional  areas  reset: 
Top  =  0. 

Right  =  0. 

Bottom  =  0. 
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The  calculation  of  the  fractional  areas  is  thus  re 
duced  to  considering  four  directions  of  an  interface  with 
respect  to  a  specified  cell  K: 

1.  entering  cell  K,  leaving  cell  KR 

2.  entering  cell  KR,  leaving  cell  K 

3.  entering  cell  K,  leaving  cell  KA 

4.  entering  cell  KA,  leaving  cell  K. 


4.4  COMPUTING  THE  FLUX  TERMS  FOR  INTERFACE  CELLS 

The  computation  of  fractional  cell  boundary  areas  is 
performed  in  order  to'  use  the  material  interface  positions 
to  define  the  mass  flux  of  materials  across  interface  cell 
boundaries.  For  pure  cells  the  mass  flux  between  two  cells 
is  ,  , 

Am  =  P.J.  A  At 

where  is  the  donor  cell  density,  is  an  average  of 

the  cell-centered  velocities  of  the  two  cells,  and  A  is  the 
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area  of  the  cell  boundary  for  which  the  mass  flux  is  being 
computed. 

For  cells  crossed  by  an  interface  this  equation  becomes 
AmN  =  V.J,  FAj^  At/S 

where  is  usually  the  donor  material  N!  density,  FA 

is  a  fractional  cell  boundary  area  for  material  N,  and  the 
time  step,  At,  is  divided  by  S,  the  number  times  the 
INFACE  subroutine  is  subcycled.  Thus  an  interface  cell  has 
mass  flux  terms  for  its  right  and  top  boundaries  associated 


with  each  material  package  in  the  grid. 


For  cells  containing  the  void  interface,  "free  surface"  cells, 
the  code  applies  a  special  rule  to  define  the  transport  den¬ 
sity,  P-pjj.  If  two  neighbor  cells  are  both  free  surface  cells, 
the  transport  density  is  chosen  in  the  usual  manner  according 
to  the  direction  of  the  flow,  i.e.,  donor  material  density. 

If  only  one  of  the  two  cells  is  a  free  surface  cell,  the 
density  of  the  material  in  the  non-free  surface  cell  is  used 
regardless  of  the  direction  of  the  flow. 
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In  the  case  drawn  above,  assume  there  are  three  material 
packages  in  the  grid.  There  are  five  non-zero  mass  flux 
terms  associated  with  cell  K.  Using  the  FORTRAN  variables 
and  assuming  M  is  the, location  in  the  mixed  cell  arrays 
associated  with  cell  K,  the  right  boundary  fluxes  would  be 


jSAMMPCl  ,M)  1 

'V  Pti 

U,p  .  2921 

■  y'j-i) 

1  2tt 

^i 

At/ s 

|SAMMP(2,M) 1 

PT2 

.375( 

■  >"^-1) 

1  2tt 

x . 

1 

At/s 

1SAMMP(3,M) 1 

Pt3 

U,j.  .3331 

2tt 

^i 

At/s 

and  the  top  boundary  fluxes  would  be 
|SAMPY(l,M)i  .417(x| 

|SAMPY(2,M)  I  0-  .  583(x? 

|SAMPYC3,M3  1=0. 

The  sign  of  these  fluxes  indicates  the  direction  of  the  flow. 
The  fractional  top  cell  boundary  area  for  the  material  of 
package  three  is  zero.  Therefore,  the  material  of  package 
three  is  not  transported  across  the  top  boundary  of  this  cell. 

4.5  CREATING  MIXED  CELLS 

4.5.1  Assigning  Storage  Location 

Initially  those  cells  containing  a  material  interface 
are  created  in  SETUP  when  the  problem  is  generated.  Later, 
as  the  interfaces  move  across  the  grid,  new  mixed  cells  are 
created  in  NEWMIX  (called  by  INFACE) .  Whenever  a  cell  boundary 
is  straddled  by  a  pair  of  consecutive  tracers,  INFACE  checks 
that  both  cells  (K  and  K")  on  either  side  of  the  boundary 
are  mixed.  If  one  or  both  are  pure  (i.e.,  MFLAG(K)  <  100, 


-  IT  At/s 

-  X?  At/s 

1-1^ 


51 


3SR-350 


or  MFLAGCK")  <  100),  subroutine  NEWMIX  is  called  and  an 

’ft 

unused  storage  location  in  the  mixed  cell  arrays  is  assigned 
to  that  cell.  (The  maximum  number  o£  mixed  cells  at  any  one 
time  is  controlled  by  the  input  variable  NMXCLS  which  should 
correspond  to  the  dimensions  of  the  mixed  cell  arrays.  If 
the  user  tries  to  generate  more  than  NMXCLS  mixed  cells, 

NEWMIX  calls  for  an  error  exit.)  Before  a  cell  K  becomes 
mixed,  its  flag,  MFLAGCK),  indicates  what  material  package 
it  belongs  to.  NEWMIX  uses  that  information,  transferred 
through  a  variable  MO  in  blank  common,  to  define  the  mixed 
cell  variables  for  the  newly  mixed  cell,  as  illustrated  by 
the  following  FORTRAN  statements. 

MO  =  MFLAGCK)  (when  pure) 

(M  =  storage  location  for  cell  K 
in  mixed  cell  arrays) 

MFLAG(K)  =  M+lOO  (flag) 

XMASS(MO,M)  =  AMX(K)  (mass) 

SIE(M0,M)  =  AIX(K)  (specific  internal  energy) 

RHO(MO,M)  =  AMX(K)/[TAU(I)*DY(J)]  (density) 

A  cell  that  becomes  "mixed”  in  INFACE  has  only  one 
material  until  subroutine  PH2  is  executed  and  the  other 
materials  are  transported  into  it.  Note  also  that  PH2  will 
not  transport  material  of  another  package  into  a  pure  cell. 

The  interface  of  the  other  package  must  first  cross  the 
cell's  boundary  in  INFACE  before  the  cell  can  receive  the 
material  of  the  other  package  in  PH2 . 


A  storage  location  M  is  not  in  use  if  RHO(l,M)  =  -1. 
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4.5.2  Defining  Free  Surface  Cells 

A  cell  containing  a  free  surface  often  contains  only 
one  material  and  yet  is  still  treated  as  a  mixed  cell  since,  . 
an  interface  crosses  its  boundaries  and  the  fluxes  across 
its  boundaries  are  a  function  of  the  fractional  areas  computed 
in  INFACE.  Of  course  it  is  possible  for, a  cell  containing 
a  void  (i.e.,  a  free  surface)  to  contain  more  than  one 
material.  A  "free  surface”  cell  K  is  flagged  by  setting 
RHO(NVOID,M)  =  1.,  where  NVOID  is  one  more  than  the  number 
of  material  packages  in  the  problem  and  M  =  MFLAG (K) -100 . 

4.5.3  Accounting  for  Subcycles 

* 

INFACE  usually  is  subcycled,  i.e.,  the  routine  is 
executed  several  times  each  cycle  using  a  fraction  of  the 
time  step  on  each  subcycle.  It  is  therefore  possible  to 
create  a  mixed  cell  after  one  or  more  subcycles  of  INFACE 
have  been  completed.  In  that  case  -the  flux  terms  (SAMMP, 
SAMPY)  for  that  cell  have  not  been  accumulated  for  the  com¬ 
pleted  subcycles  and  need  to  be  computed  before  adding  in 
the  flux  terms  for  the  current  and  subsequent  subcycles. 
Therefore,  when  INFACE  has  completed  at  least  one  subcycle, 
NEWMIX  calls  FLUX  after  setting  up  the  storage  for  the  new 
mixed  cell. 

4.6  CREATING  PURE  CELLS 

4.6.1  Defining  Material  that  Remains  in  the  Cell 

After  the  fractional  areas  associated  with  each 
material  package  have  been  computed,  INFACE  searches  for 
cells  that  were  interface  cells  on  the  previous  cycle  but 

•k 

Subcycling  INFACE  minimizes  transport  noise. 
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no  longer  are  cut  by  an  interface.  Such  a  cell  has  become 
pure  and  its  flag,  MFLAG,  is  made  negative  to  signify 
that  fact  until  the  transport  has  been  completed  in  PH2.  To 
determine  what  material  package  a  newly  pure  cell  K  belongs 
to,  INFACE  first  looks  for  a  neighbor  that  is!  pure  and 
assumes  cell  K  will  belong  to  the  same  package  as  one  of 
its  pure  neighbors,  e.g.  in  the  following  diagram, 


the  cell  below  is  pure  CMFLAG(KB)  =  2),  therefore  cell  K 
has  become  a  pure  cell  belonging  to  package  2.  In  case  all 
four  neighbors  are  mixed,  the  cells  KB  and  ' KL  are  examined 
as  illustrated  in  the  following  example. 
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Since  cell  K  is  not  crossed  by  an  interface  it  follows  that 
if  the  fractional  area  of  the  top  of  cell  KB  for  package 
2  is  the  total  cell  boundary  area  then  cell  K  must  belong 
to  package  2.  The  cell  on  the  left  of  cell  K  is  used 
analogously  if  cell  K  is  in  the  bottom  row  of  the  grid. 


4.6.2  Evacuating  Materials 


To  signify  a  material  interface  has  passed  out  of  cell  K, 
the  density  of  that  material  is  set  to  zero.  For  example,  con¬ 
sider  the  figure  above  assuming  there  are  three  material 
packages  in  the  grid;  the  densities  of  cell  K  are  redefined 
until  the  end  of  the  cycle  as  follows: 


M  =  MFLAGCK)  -  100 
RH0C1,M)  =0. 
RHO(2,M)  f  0. 
RHO(3,M)  =  0. 
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These  zero  densities  are  later  used  as  a  signal  to  adjust  the 
fluxes  of  materials  1  and  3  so  they  will  be  exactly  evauated 
from  cell  K.  This  adjustment  of  the  fluxes; usually  occurs 
twice  in  INFACE,  once  during  the  first  subcycle  (for  cells 
that  become  pure  during  the  first  subcycle)  and  once  after 
all  subcycles  are  completed  (for  cells  that  become  pure  after 
the  first  subcycle).  For  the  first  case  the: flux  terms  for 
mixed  cells  (SAMMP,  SAMPY)  are  saved  from  the  previous 
cycles  to  indicate  the  direction  of  the  flow  and  the  direction 


of  the  evacuation.  In  the  following  example;  t  repre¬ 
sents  time,  and  ML  and  MB  are  locations  associated  with  the 


mixed  cells  on  the  left  of  and  below  cell 


K, 


respectively: 


t  flux  terms 
1 

SAMMP (1,M)  >  0. 
SAMPY (I, M)  >  0. 
SAMMP (1, ML)  =  0. 
SAMPY (1, MB)  =  0. 


At  t  ,  the  t  flux  terms  indicate  that  any  mass  of 
2  1 

package  1  still  in  cell  K  should  be  evacuated  out  the 
right  and  top  boundaries,  not  the  left  or  bottom.  These 
evacuation  procedures  are  used  for  a  cell  that  has  lost 
one  interface  but  is  still  cut  by  another  interface  (so 
is  still  mixed)  as  well  as  for  cells  that  have  become 
pure.  I 
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4.7  ADVANCING  TRACER  PARTICLES 

Each  material  is  circumscribed  by  a  series  o£  massless 

tracer,  particles,  which  are  propagated,  each  time  step  a  dis- 
* 

tance 

Ax .  =  u .  At 
1  1 

where  is  a  local  average  velocity  vector  for  the 

continuum,  determined  by  an  area  overlap  method  which  gives 
weight  to  velocities  in  the  surrounding  cells.  Specifically, 
a  rectangle  of  cell  dimensions  is  superimposed  on  the  particle 
to  be  moved  and  then 


Since  the  tracer  coordinates  are  in  grid  line  units,  this 
distance  is  converted  from  centimeters  to  cell  units. 
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where  wfL)  is  the  cell-centered  velocity  of  the  overlapped 
cell  L  and  Aj^  is  the  overlap  area.  This  procedure  gives 
a  spatially  continuous  velocity  field  for  particle  propagation. 

The  cell-centered  tracers  (XP,  YP) ,  used  only  for 
plotting  the  movement  of  mass  within  the  package  boundaries, 
are  moved  by  the  same  method. 


These  cell-centered  velocities  have  been  updated  by  PHI 
(pressure  effects)  and  PH3  (strength  effects)^  for  this 
time  step. 
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V.  THE  PRESSURE  AND  DENSITY  ITERATION 
FOR  MIXED  CELLS 


5.1  GENERAL  METHOD 

For  ordinary  one-material  cells,  it  is  possible  to 
compute  the  pressure  directly  from  the  material  density  and 
specific  internal  energy.  For  cells  containing  more  than 
one  material,  it  is  necessary  to  make  additional  assumptions 
or  observations  in  order  to  determine  cell  pressure.  In  the 
present  method  we  determine  the  cell  pressure  for  mixed  cells 
by  an  iteration  procedure.  Specifically,  the  densities  of 
the  various  materials  within  the  cell  are  varied,  subject 
to  the  condition  that  the  cell  be  exactly  filled  by  the 
masses  therein,  until  the  individual  pressures  converge  to 
a  common  value,  taken  to  be  the  cell  pressure.  Specific 
internal  energies  are  held  constant  during  the  iteration. 

This  process  gives,  in  addition  to  the  cell  pressure,  the 
densities  of  the  individual  materials  for  subsequent  use 
in  the  Phase  1  and  Phase  3  energy  partition  and  in  the 
Phase  2  transport  calculation. 

The  procedure  for  determining  the  pressure  and  the 
densities  in  a  mixed  cell  is  actually  divided  into  two  parts, 
a  pre- iteration  calculation  to  assure  that  the  cell  is  filled 
exactly,  and  an  iteration  calculation  to  converge  on  the 
desired  P's  and  p's. 

5.2  PRE- ITERATION  CALCULATION 

In  general,  the  masses  and  densities  which  exist  at 
the  beginning  of  the  calculation  do  not  exactly  fill  the  cell, 
since  cell  masses  were  changed  in  the  Phase  2  calculation  of 
the  previous  cycle.  The  pre-iteration  calculation  fills  the 
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cell  exactly,  taking  account  o£  the  (computed)  compressi¬ 
bilities  of  the  various  materials  in  the  necessary  alteration 
of  densities.  Stepwise,  the  code 

1.  Computes  from  the  old  value  of 

for  each  material  within  the  cell  ;(E^  is 
assumed  constant  throughout) . 

2.  Varies  the  p-  by  one  percent  to  compute 

1  ■  ■  j  ; 

new  P . .  .  ■ 

1  .  ;  ■ 

3.  Determines  the  constant-energy  compressi¬ 
bilities,  called  C?,  for  each  material 

C?  =  AP./Ap. 

1  i'  ■  1 

where  the  AP^^  and  Ap^  are  given  from 
the  results  of  (1)  and  (2) . 

4.  Normalizes  (see  proof  below)  to  fill  cell 

exactly:  j 


where 


and 


AV^  =  -AP/h? 

hi  =  (pj  c.)" 


AP  =  (VOL-VCELL) 

2 
i 


VOL 
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=  mass  of  constituent  i. 
VCELL  =  volume  of  the  cell. 


With  the  above  increments  in  specific  volume, 
then,  the  new  specific  volumes  are  given,  by 


V. 

.1 


i,old 


+  AV. 
1 


then. 


Pi  ■  l/Vi 


gives  the  desired  new  densities,  causing  the 
cell  to  be  filled  exactly. 


It  is  easy  to  see  that  the  calculation  defined  above 
causes  the  cell  to  be  filled  exactly.  To  show  this. 


m 


“l/Pi.old  -  "'l 


m .  /o .  , j 

1  , old 


m^CVOL-VCELL) 


m . 


h? 


Summing , 


2jm.V.  =  VOL 
■  1  1 
1 


-  (VOL -VCELL) 


h^ 

1 


m . 
2 

h. 


=  VCELL 
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showing  that  the  new  volumes  of  the  constituents  add 

up  to  the  total  cell  volume,  as  desired,  ' 

I 

5.3  ITERATION  CALCULATION 

-  I 

Given  E-)  for  materials  i  =  1,M  find 

cell  pressure  F=P  =P  =...=Pwby  varying  the  V- 

12  M  1 

until  the  P^  are  equal  within  some  specified  accuracy, 

subject  to  the  conditions  that  the  cell  remain  exactly 

filled  and  that  the  E.  are  constants. 

1 

The  equation  for  volume  conservation  is 


^  m^  AV^  =  0 

and  the  equation  that  causes  material  i  to  Undergo  a  change 
in  specific  volume  AV^,  such  that  its  pressure  change  from 
its  current  value  P^  =  to  a  value;  P^^^,  common 

to  all  materials  in  the  cell  at  the  end  of  the  iteration  cycle, 
is 


These  equations,  for  i  =  1  ...  M  can  be  regarded  as  M+1 

~N+1 

equations  for  the  M+1  quantities  AV^  and  P  .  The 
other  quantities  in  the  equations  are  either  known  constants 
(m^)  or  are  updated  each  cycle  of  the  iteration  [the  P^  and 
the  (9V./3P.)p  ]  and  are  taken  to  be  constaiits  while  solving 

these  equations  for  AV-  and  P^  .  The  solutions  are 


^+1 


Si: 
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where 


h.  =  (3Pi/SVj) 


That  these  equations  are  solutions  to  the  given  equations  can 
be  verified  by  direct  substitution. 

The  iteration  would  be  exact  (single  cycle  convergence^ 
if  the  input  coefficients  were  constants,  since 

the  solution  does  not  involve  additional  approximations.  In 
the  iteration,  these  coefficients  are  refined  each  step  of 
the  iteration  by  making,  use  of  the  P^,  points  which  were 
generated  during  the  previous  step.  Further,  some  pains 
are  taken  to  determine  realistic  starting  values  of  these 
coefficients,  in  the  pre- iteration  calculation,  in  the  belief 
that  this  procedure  will  hasten  convergence  sufficiently  to 
reduce  total  computing  time. 

Stepwise,  the  iteration  calculation  is  as  follows: 

1.  Compute  new  and  AVi  by  the  above 

formulas,  using  known  values  of  and 

1 

2.  Compute  new  by  adding  the  above  AV^ 

to  the  old  specific  volumes. 

3.  Compute  new  P^^^,  using  the  updated  . 

■  4.  Compute  refined  values  of  the  coefficients 
^9V^/9P^)g  by  using  this  last  point  [from 
(2)  and  (s)  above]  and  that  previous 
point  (two  such  points  are  saved  for  each 
material  during  the  iteration)  which  is 
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I 

closest  in  pressure  to  this  last  point, 

i.  •  6  •  I  1 


where  P^,  is  the  closest  previous  point. 

5.  Return  to  (1) »  replacing  the  old  values  of  the 
and  the  with  the  iupdated 

values.  ^  ; 

The  iteration  is  completed  when  the  P^  all  equal  P  to 
some  preset  accuracy.  In  most  experience  to  date,  convergence 
has  been  obtained  in  two  or  three  steps,  where  a  convergence 
criterion  |P^  -  P|  <  10"^  P  was  used.  i 

5.4  EQUATION  OF  STATE  MODIFICATIONS  FOR 

THE  ITERATION  PROCEDURE  ' 

The  Tillotson  Equation  of  State  Formulation  described 
in  Section  2.3,  has  been  used  in  most  work  to,  date  with  HELP 
and  is  sufficiently  general  to  provide  a  satisfactory  descrip¬ 
tion  of  most  media  which  is  initially  condensed  and  which  may 
vaporize  as  a  result  of  heating. 

For  pure  cells,  certain  physically  meaningless  com¬ 
puted  pressures  are  excluded  after  they  are  calculated,  i.e., 
negative  pressures  corresponding  to  expansions  to  densities 
less  than  some  input  number.  Also,  sufficiently  expanded 
states,  for  some  materials,  can  lead  to  the  abnormal 
(dP/dV  >  0)  slope  of  a  constant  energy  curve  with  E  <  E^. 
These  meaningless  states  are  avoided  by  requiring,  in  the 
pure  cell  formulation,  that  the  quantity  ^Ay  +  By^j  ,  where 
y  =  (p/p|j  -  l),  is  replaced  by  its  mininum  value  (a  negative 
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pressure)  if  p  is  less  than  the  density  corresponding  to 
this  minimum  value. 

For  mixed  cells,  to  assure  convergence  of  the  present 
iteration,  it  is  desirable  to  have  a  continuous  equation  of 
state  and  for  9P/3V  to  be  everywhere  negative.  To  this  end 
the  term 


Ay  +  By  ^ 


is  replaced  by 


Ay  +  hCV)By^ 


where 

1.  hCV)  =  +1  for  V  <  V  compressed 

~  0 

states , 

2.  h(V)  decreases  linearly  to  -1  in 
the  interval 


V 

0 


<  V  <  V 

0 


2B+A 

TFX 


and 


3.  hCV) 


-1 


for 


V  >  V 
—  0 


2B+A 


The  effect  of  these  modifications  is  to  change  B  to  -B 
in  a  continuous  manner  and  hence  insure  that  3P/SV  is 
always  negative.  Negative  pressures  that  may  result  after 
the  iteration  is  completed  are  set  to  zero,  since  tensile 
states  in  a  mixed  cell  are  assumed  not  to  exist.  In  view 
of  this  latter  assumption,  the  indicated  equation  of  state 
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changes  do  not  normally  affect  the  physical  result,  since 
they  alter  only  states  that  would  ordinarily  be  tensile. 

The  purpose  of  the  changes  is  to  make  the  equation  of  state 
well-behaved,  so  that  the  iteration  converges  properly. 
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VI.  SAMPLE  CALCULATIONS 


6.1  INTRODUCTION 


HELP  was  written  as  part  of  a  theoretical  study  of 
the  armor  penetration  process,  but  is  also  applicable  to 
a  wide  range  of  problems  in  continuum  mechanics.  The  code, 
sometimes  with  modifications  or  additional  routines,  is  now 
used  to  treat  problems  in  gas  dynamics,  underground  explosions, 
explosive-metal  interactions,  shaped  charges,  radiation 
hydrodynamics  and  a  variety  of  impact  applications.  The  pro¬ 
gram  is  especially  well  suited  for  solving  problems  involving 
extreme  material  distortions,  in  which  there  may  be  more  than 
one  material  in  the  flow  field.  In  the  following  pages, 
selected  results  are  given  from  several  problems  solved  with 
HELP. 
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6 . 2  EXPANSION  OF  SPHERICAL  SHELLS 

The  purpose  o£  this  test  problem  was  to  affirm  that 
the  code  would  preserve  the  spherical  symmetry  during  an 
explosion  of  a  hot  spherical  gas  bubble  intO' a  surrounding 
medium  of  concentric  shells.  The  test  is  meaningful  because 
the  calculation  is  performed  as  a  2-D  problem  in  axisymmetric 
flow,  so  that  any  difference  in  treating  the  axial  and  radial 
motions  should  appear  as  an  asymmetry. 

The  initial  conditions  are  those  indicated  in  the 

Y  E 

0 

1.5  10^® 

1.4  0. 

1.2  0. 

1.3  0. 

;s) 


The  four  times,  corresponding  to  automatic  plots  of  the 
tracer  positions  in  Fig.  3,  are  0.,  0,05,  0.13,  and  0.23  psec. 

As  seen  from  these  plots,  the  solution  retains  the 
desired  spherical  symmetry  to  the  latest  times  computed. 
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6.3  IRON  PROJECTILE  IMPACTING  AN  ALUMINUM  TARGET 

This  problem,  discussed  in  detail  in  Ref.  6,  is  the 
impact  of  an  iron  projectile  on  an  aluminum  target.  The 
impact  velocity  is  9  ^  10**  cm/sec  and  the  shear  yield  strengths 
of  the  iron  and  aluminum,  respectively ,. are  6!  kilobars  and 
3  kilobars. 

1  cm 

1.15  cm 


Times  presented  in  the  plots  of  Fig.  4  are  0,2.5,  5.0  and 
10.0  ysec.  In  the  last  frame,  the  10  ysec  plot  is  repeated 
with  the  positions  of  the  interior  (passive)  tracers  shown 
instead  of  the  Eulerian  grid  lines.  In  all  plots  the  out¬ 
lines  of  the  target  and  projectile  are  obtained  in  the  plot 
routine  by  drawing  straight  lines  between  the 'free  surface 
and  interface  tracers. 
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Fig.  4--Iron  projectile  impacting  an  aluminum  target. 
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6 . 4  IMPACT  OF  A  CONE  ONTO  WATER 

To  simulate  an  experimental  study  of  impact  into  water, 
a  conical  tungsten  projectile  (solid  near  tip,  thin-walled 
behind)  impacts  into  water  at  2.4  x  10®  cm/se‘c.  The  yield 
strength  of  the  tungsten  is  10.5  kilobars.  Zoning  is  variable 
(0.2  cm  zones  near  impact  center,  up  to  0.5  cm  in  far  field). 

Automatic  plots  of  the  configuration  are  given  at  six 
times  in  Fig.  5.  Times  are  0,  2.5,  11,  18,  2'3  and  30  ysec. 

The  solid  tip  is  initially  2.2  cm  in  length.  As  seen  from 
these  plots,  the  tip  is  substantially  eroded  by  the  penetra¬ 
tion  process.  At  the  latest  times  in  the  calculation  the 
water  splash  is  beginning  to  interact  with  the  thin  tungsten 
walls. 


1 
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6.5  METEOROID  IMPACT  INTO  LUNAR  SURFACE 

A  stony  meteoroid,  10  km  in  diameter,  impacts  the 

moon  at  3  X  10®  cm/sec.  The  lunar  surface  is  taken  to  be 

i 

basalt  with  20  percent  porosity. 

Five  times,  0,  0.4,  0.6,  0.9  and  2.5  seconds  are 
shown  in  Fig.  6.  At  the  second  time,  0.4  sec,  the  shock  wave 
has  just  reached  the  upper  surface  of  the  meteoroid.  In  subse¬ 
quent  plots  the  upper  surface  is  seen  to  fold; back  upon  the 
axis  of  symmetry,  the  meteoroid  itself  (see  0.9  sec)  forming 
a  thin  layer  at  the  bottom  of  the  growing  crater.  At  the 
later  times  shown,  the  shock  wave  into  the  moon  is  apparent 
as  a  region  in  which  the  tracers  are  relatively  dense. 

Detailed  analysis  of  the  output  gives  also  the  angular 
and  velocity  distribution  of  the  debris  and  the  regions  near 
the  center  of  impact  in  which  the  lunar  material  is  vaporized 
or  melted  by  the  shock  heating. 
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Fig. 


6--Meteoroid  impact  into 


luna_ 
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6.6  JET  FROM  A  COLLAPSING  METAL  LINER 


A  copper  liner,  a  shortened  version  o£  that  in  the 
105-mm  shaped  charge,  is  imploded  on  axis.  Initial  velo¬ 
cities  are  along  the  liner  wall  at  3.3  x  10^  , cm/sec  except 
near  the  axis,  where  the  velocity  drops  linearly  to  zero 
over  the  last  0.3  cm.  The  initial  thickness  to£  the  liner 

I 

wall  is  0.27  cm  and  the  Eulerian  grid  consisted  of  0.05  cm 
square  cells. 

i 

Times  shown  in  the  plots  are  0,  0.7,  1.5  and  2.6  ysec 
At  the  later  times,  the  copper  jet  is  seen  to  be  well  formed 
with  the  main  part  o£  the  impinging  streams  coalescing  into 
the  slug.  This  test  problem  is  part  of  a  general  study 
being  conducted  by  M.  L.  Gittings,  R.  T.  Sedgwick  and  the 
present  authors,  for  predicting  the  performaiice  of  shaped 
charges. 
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6.7  ARTILLERY  SHELL  EXPLODING  IN  AIR  I 

i 

An  in-flight  10 5 -mm -diameter  artillery  |shell,  3000  £t/sec, 
is  exploded  in  air.  The  shell  is  filled  with  Composition  B 
explosive  and  has  a  0.5-cm  steel  wall  and  a  sodid  steel  nose 
section. 

Explosive  detonation  is  assumed  to  begin  at  the  front 
and  proceed  to  the  base  of  the  projectile.  Various  times, 

0,  30,  47,  74,  124  and  203  ysec  are  shown. in  the  plots  of 
Fig.  8,  where  the  steel-air  interfaces  are  plotted  as 
asterisks  and  the  explosive  and  air  tracers  arp  represented 
by  periods.  A  rigid  wall  is  present  at  the  boittom  boundary  of 
the  plots. 

! 

The  air  blast  wave  is  seen  to  have  a  smooth  and  nearly 
hermispherical  shape  at  intermediate  times.  At  the  203  ysec 
time,  however,  the  steel  fragments  in  the  forward  half  of  the 
flow  have  nearly  overtaken  the  air  shock  wave.  At  later  times 

still,  these  fragments  penetrate  the  main  blast  wave,  which  is 

I 

decelerating  due  to  divergence.  A  more  detailed  description 

[71 

of  this  test  problem  is  given  in  a  recent  publication.  ■' 


■I 
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VII.  GUIDE  TO  THE  FORTRAN  PROGRAM 

7.1  INTRODUCTORY  REMARKS 

The  HELP  code  is  basically  the  RPM  code  with  the  added 
capability  of  treating  material  interfaces  and  of  giving  an 
elastic-plastic  description  of  material  strength.  As  in  RPM 
the  equations  of  motion  are  integrated  in  three  phases,  PHI, 
PH2  and  PH3,  which  account  for  the  effect  of  pressure,  trans¬ 
port  and  shear  stresses  respectively.  Subroutine  INFACE  has 
been  added  to  move  the  material  interfaces  and  to  compute  the 
mass  flux  of  cells  cut  by  one  or  more  interfaces.  An  itera¬ 
tion  to  equilibrate  pressures  in  mixed  cells  has  been  added 
to  CDT,  and  EQST  gives  pressure  as  a  function  of  density 
and  energy  for  nineteen  materials.  The  GLUE  routine  treats 
partially  filled  free  surface  cells  and  cells  that  have  been 
over  accelerated  by  the  Phase  1  effects  of  pressure.  These 
and  the  more  peripheral  routines  of  HELP  are  described  in  more 
detail  in  the  paragraphs  below.  A  general  flow  diagram  of 
HELP  is  given  in  Fig.  9,  and  references  to  subroutines  are 
summarized  in  Table  2. 

7.1.1  ADDTCR 

Material  tracers  can  be  added  in  a  specified  region  at 
any  point  in  the  calculation  by  defining  input  flags  (NADD, 
MINX,  MINY,  MAXX,  MAXY) .  When  NADD  >  0,  ADDTCR  is  called  from 
EDIT  on  the  first  cycle  of  a  run.  ADDTCR  adds  at  most  the 
specified  number  (NADD)  of  tracers  between  each  existing  pair 
of  tracers  in  the  region  but  „does  not  create  a  greater  density 
of  tracers  than  that  initially  specified  by  NTRACR. 
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For  pure  hydro 
calculations 
PH3  is  omitted* 


Computes  effects 
of  deviatoric  and. 
hoop  stresses  to 
update  velocities 
and  energies  . 


Computes  fluxes 

FLUX 

of  materials  in 

NEWMIX 

mixed  cells. 

NEWRHO 

Advances  tracer 
particles. 

Ends  Cycle 


Accounts  for  mass 
flux  and  asso¬ 
ciated  transports 
of  momentum  and 
energy . 


Fig.  9--General  flow  diagram  of  HELP  program. 
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TABLE  2 

SUBROUTINE  REFERENCES 


Name 

Called  From 

Calls 

ADDTCR 

EDIT 

CARDS 

INPUT 

CDT 

HELP,  EDIT 

ERROR,  EQST 

COMPRS 

REZONE 

EDIT 

HELP,  ERROR 

ADDTCR,  CDT,  ERROR,  MAP, 
REZONE 

EQST 

CDT 

ERROR 

CDT,  EDIT,  INFACE, 
INPUT,  NEWMIX,  NEWRHO, 
PH2 

EDIT 

FLUX 

INFACE,  NEWMIX 

GLUE 

PHI 

HELP 

CDT,  EDIT,  INFACE,  INPUT, 
PHI,  PH2,  PH3 

INFACE 

HELP 

ERROR,  FLUX,  NEWMIX, 
NEWRHO 

INPUT 

HELP 

CARDS,  ERROR,  SETUP 

LOCIJ 

SETUP 

-- 

MAP 

EDIT 

-  - 

NEIVMIX 

INFACE,  REZONE  , 

ERROR,  FLUX 

NEWRHO 

INFACE 

ERROR 

PHI 

HELP 

GLUE 

PH2 

HELP 

ERROR 

PH3 

HELP 

STRNG 

PROPRT 

SETUP 

RE ZONE 

EDIT 

COMPRS,  NEWMIX 

SETUP  , 

INPUT 

LOCIJ,  PROPRT 

STRNG 

PH3  ■ 

— 
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7.1.2  BLOCK  i 

A  block  data  element  that  defines  normal  density  and 

t 

the  bulk  sound  speed  of  nineteen  materials.  | 

I 

7.1.3  CARDS 

— I— j 

Most  of  the  input  parameters  stored  in  blank  common 
are  read  by  subroutine  CARDS.  The  format  of  these  input  cards 
is  described  in  Section  7.2.2. 

7.1.4  COT  , 

A  principal  function  of  this  routine  is  to  compute  a 
time  step  which  ensures  stability  of  the  finite  difference 
equations.  This  is  done  by  finding  the  cell  with  the  minimum 
ratio,  D/o).  Here  D  denotes  the  minimum  of  the  cell's  dimen¬ 
sions  (cm),  and  co  denotes  the  sum  of  the  cell's  maximum 
velocity  component  and  its  sound  speed  (cm/sec) .  For  a 
polytropic  gas  the  sound  speed  is  computed  as  ’  /yp/p  and 
for  other  materials  by  the  approximate  relation  C  =  C^  +  1” 
where  P  is  the  pressure  in  the  cell.  The  coefficient  S’ 
is  an  input  parameter  and  is  described  more  fully  in  Section  7.3 

C  is  defined  in  BLOCK  for  nineteen  materials  as  /A/p  , 
a  0 

where  A  is  a  coefficient  in  the  solid  equation  of  state. 

In  a  mixed  cell  the  sound  speed  is  given  by  a  ;mass  weighted 
average  of  the  sound  speeds  of  the  materials  in  the  cell. 

Each  cycle  CDT  prints  the  column  and  row  (I,J)  of  the  cell 
controlling  the  time  as  well  as  the  maximum  sound  speed  plus 
velocity  in  the  grid  (MAXCUV) ,  the  maximum  velocity  in  the 
grid  (MAXUV)  and  the  velocity  and  pressure  cutoff  values. 

Another  function  of  CDT  is  to  equilibrate  the  pressures 

1 

of  materials  in  mixed  cells  using  an  iteration  scheme  that 
adjusts  the  material  densities.  A  detailed  discussion  of  this 
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iteration  method  is  found  in  Section  V  .  The  pressures 
for  pure  cells  are  also  updated  in  CDT  by  a  call  to  the 
equation  of  state  subroutine,  EQST,  with  the  density,  energy 
and  material  code  number  of  the  cell  CRHOW,  ENERGY,  N)  passed 
through  blank  common. 

7.1.5  COMPRS 

{ 

COMPRS  is  called  by  subroutine  REZONE  to  create  a  new 
cell  by  combining  two  adjacent  cells  into  one  and  defining 
the  velocity  and  energy  of  the  new  cell  so  that  energy  and 
momentum  are  conserved.  If  both  cells  are  pure  then  the  new 
cell  will  also  be  pure.  If  one  or  both  are  mixed,  the  new 
cell  will  be  mixed. 

7.1.6  EDIT 


The  periodic  printing  and  writing  on  the  restart  tape 
are  executed  by  subroutine  EDIT.  The  frequency  of  printing 
and  tape  dumps  is  controlled  by  input  parameters  described  in 
Section  7.2.3. 

On  every  print  cycle  EDIT  prints  the  mass,  total  energy, 
internal  energy,  kinetic  energy,  axial  and  radial  momentum, 
and  plastic  work  for  each  material  package  as  well  as  for  the 
entire  grid.  The  changes  in  energy  due  to  evaporation  and 
losses  out  boundaries  are  also  accounted  for.  The  coordinates 
of  the  material  tracers  circumscribing  each  package  are  printed 
in  cell  units.  Summary  graphs  of  the  compression,  pressure, 
velocities  and  internal  energies  are  printed  if  the  input 
parameter,  MAPS,  is  non-zero.  And  finally,  the  pressure,  velo¬ 
cities,  internal  energy,  compression  and  stress  deviators  for 
all  cells  in  the  active  grid  are  displayed  on  "long"  EDIT  prints 
and  for  each  cell  in  the  axis  column  on  "short"  EDIT  print. 
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Another  important  function  of  EDIT  is  to  compute  the 
relative  error  in  the  total  energy  of  the  grid  when  summed 
over  all  the  cells,  and  to  check  that  this  error  does  not 
exceed  a  limit  specified  by  the  input  variable  DMIN, 

EDIT  also  senses  when  execution  shoul<i  stop  and  sets 
the  exit  flag,  WFLAGL. 

7.1.7  EQST 

The  equation  of  state  constants  for  nineteen  materials 
are  stored  in  DATA  statements  in  EQST.  (These  constants  are 
displayed  in  Table  1.)  Given  a  code  number  for  specifying 
the  material,  a  density  and  an  energy,  EQST  computes  a 
pressure.  A  more  detailed  discussion  of  the  equation  of 
state  is  given  in  Section  2.3.  EQST  is  called  by  CDT  to 
compute  pressures  of  pure  cells  as  well  as  of  materials  in 
mixed  cells.  ' 

7.1.8  ERROR  1 

This  subroutine  is  called  in  the  case  of  certain  error 
conditions  tested  on  by  the  code,  e.g.,  the  energy  sum  error 
exceeding  the  specified  limit.  ERROR  prints  a  message  identi¬ 
fying  the  general  location  of  the  error  condition,  lists  the 
Z-block  variables  (first  150  words  of  blank  common) ,  calls 
EDIT  to  do  a  long  print  and  a  tape  dump,  and  then  calls  exit. 

7.1.9  FLUX  ! 

The  mass  flux  across  the  right  and  top  boundaries  of 
a  mixed  cell  for  each  material  in  that  cell  is  computed  in 
FLUX.  These  fluxes  are  stored  in  the  SAMPY  and  SAMMP  arrays 
and  executed  in  subroutine  PH2  which  does  the  actual  transport 
of  material  across  all  cell  boundaries.  FLUX  uses  the  position 
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of  the  material  interfaces  to  compute  the  fluxes  of  the 
various  materials  in  the  cell.  See  Section  4.2 
for  a  more  complete  description  of  transport  from  mixed 
cells . 

7.1.10  GLUE 

Free  surface  cells,  or  cells  of  very  small  mass,  can 
be  over-accelerated  by  the  pressure  forces  (PHI)  or  the 
strength  forces  CPH3) .  See  the  discussion  in  Section  3.5. 

GLUE  is  called  by  PHI  after  the  effects  of  pressure  have 
been  computed  and  after  PHI  searches  for  cells  with  un¬ 
realistic  velocities  and  makes  the  flags  (MFLAG)  of  these 
cells  negative. 

GLUE  combines  or  "glues"  a  cell  with  a  negative  flag 
or  a  free  surface  cell  to  its  adjacent  neighbor  with  the 
highest  pressure.  For  each  velocity  component,  GLUE  computes 
a  new  velocity  for  the  two  cells  which  conserves  momentum. 

The  new  specific  internal  energy  for  the  cells  is  then 
chosen  to  conserve  total  energy.  The  overall  effect  is  to 
give  the  two  cells  the  same  velocity  and  specific  internal 
energies.  The  negative  flags  are  made  positive  as  the 
cells  are  processed  by  GLUE. 

7.1.11  HELP 


The  overall  cycling  of  the  calculation  is  controlled 
by  HELP,  as  shown  in  Fig.  9.  It  is  the  main  routine  of  the 
code  and  calls  INPUT,  CUT,  EDIT,  PHI,  PH3 ,  INFACE  and  PH2 ,  in 
that  order.  When  needing  to  see  the  effects  of  each  phase 
of  the  calculational  cycle,  HELP  will  respond  to  the  input 
parameter  INTER  and  call  EDIT  on  print  cycles  after  PHI  and 
PH3  as  well  as  after  CDT.  HELP  also  calls  EXIT  on  the  normal 
cessation  of  the  calculation. 
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7.1.12  INFACE  i 

I 

The  tracer  particles  that  circumscribe  each  material 
package  are  moved  in  INFACE.  These  particles  when  connected 
by  straight  line  segments  are  used  by  INFACE  to  determine 
where  the  interfaces  cut  cell  boundaries  and  to  compute  the 
flux  of  materials  across  these  partial  cell  boundaries.  For 
a  given  cycle,  INFACE  senses  which  cells  become  mixed  and 
which  become  pure.  Section  VI  gives  a  detailed  description 
of  the  logic  and  conventions  employed  by  INFACE. 

i 

7.1.13  INPUT  . 


Instructions  for  running  problems  are  interpreted  by 
INPUT  which  can  either  start  or  restart  a  calculation.  It 
calls  SETUP  and  CARDS,  as  necessary,  to  prescribe  the  initial 
conditions  and  to  read  the  input  deck.  It  can  also  read  a 
HLPCLM  tape  for  problems  whose  initial  conditions  cannot  be 
described  by  SETUP.  - 

7.1.14  LOCIJ  ' 

This  routine  is  used  to  identify  the  column  or  row 
location  of  a  tracer  point  whose  coordinates  are  in  centimeter 
units.  SETUP  calls  LOCIJ  when  converting  the  tracer  coordinates 
to  cell  units  for  the  void  package  as  well  as  when  defining 
the  boundaries  of  rectangular  packages. 

.  i 

7.1.15 

This  subroutine  is  called  by  EDIT  when  the  input  para¬ 
meter  MAPS  is  non-zero.  It  displays  the  properties  of  each 
cell  in  the  active  grid  using  an  alphabetic  scale.  [On 
cycle  0  it  presents  the  entire  (IMAX  by  UMAX)  grid.]  One 
obtains  contour  maps  of  the  compression,  pressure,  radial  and 
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axial  velocities,  and  internal  energy.  Since  the  compression 
of  mixed  cells  is  not  computed,  asterisks  are  displayed  for 
mixed  cells  in  the  compression  map.  The  scale  of  each  map  is 
adjusted  according  to  the  current  minimum  and  maximum  values 
of  each  property. 

7.1.16  NEWMIX 

When  a  pure  cell  K  becomes  mixed,  NEWMIX  assigns  to 
it  a  location  M  in  the  mixed  cell  arrays  (SIE,  XMASS,  RHO, 
etc.)  The  flag  of  cell  K  is  redefined  so  as  to  associate 
the  K  and  the  M  indices  (MFLAGCK)  =  M  +  100).  The  fact 
that  MFLAG(K)  is  greater  than  100  indicates  cell  K  is  a 
mixed  cell.  Otherwise,  if  cell  K  were  pure,  MFLAGCK)  would 
be  equal  to  the  number  of  the  package  that  contained  cell  K. 

NEWMIX  is  called  from  INFACE  which  is  usually  subcycled. 
If  cell  K  becomes  mixed  after  the  first  subcycle  of  INFACE, 
NEWMIX  calls  FLUX  to  bring  the  definition  of  the  flux  variables 
for  cell  K  to  the  current  time.  (See  also  Section  4.5  for 
further  discussion  of  mixed  cells.) 

7.1.17  NEWRHO 

When  the  material  N  interface  first  enters  a  cell, 
INFACE  calls  subroutine  NEWRHO  to  define  the  density  of 
material  N  for  cell  K  by  looking  at  the  density  of 
material  N  in  the  neighboring  cells  closest  to  the  interface. 
(This  density  value,  stored  in  the  RHO  array,  is  updated  in 
CDT  when  pressures  are  equilibrated  at  the  beginning  of  each 
cycle . ) 

7.1.18  PHI 

The  effect  of  the  pressure  gradient  in  updating  the 
velocities  and  the  internal  energy  is  computed  here.  (The 
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numerical  method  is  described  in  detail  in  Section  2.2.2.) 
Cells  with  small  masses  that  are  over-accelerated  by  these 
pressure  effects  are  flagged,  and  subroutine’ GLUE  is  called 

to  give  more  realistic  velocities  to  these  cells. 

1 

I 

7.1.19  PH2 


Mass  transport  and  the  associated  flux  of  momentum  and 
energy  are  accounted  for  in  PH2.  (The  numerical  method  is 
described  in  detail  in  Section  2 . 2 .4 . )  After  the  transport 
has  been  completed,  PH2  corrects  for  negative  mass  quantities 
in  mixed  cells.  Before  printing  the  symbolic  map  which  dis¬ 
plays  the  material  package  numbers  and  mixed  cell  locations, 

PH2  redefines  the  flags  of  cells  that  have  become  pure. 

i 

7.1.20  pm  ]  ■ 

The  deviator  stresses  acting  on  each  cell  edge  and  the 
hoop  stress  are  determined  in  PH3  and  the  resulting  velocity 
and  internal  energy  increments  are  computed.  (Details  are 
given  in  Section  2.2.3.)  PH3  calls  STRNG  to  compute  the 
strength  of  the  material  in  a  cell.  PH3  can  be  bypassed  by 
setting  CYCPH3  =  -1.,  when  the  effects  of  strength  are  omitted. 

7.1.21  PROPRT 

Subroutine  PROPRT  is  called  by  SETUP  to  define  the  flag 
and  mixed  cell  variables  for  a  cell  on  the  edge  of  a  rectan¬ 
gular  package.  ’  I 

7.1.22  RE ZONE 

To  add  more  material  and  diminish  the  number  of  zones 
in  the  active  grid  the  REZONE  routine  can  lump  two  cells  into 
one  or  four  cells  into  one,  and  add  columns  and  rows  to  main¬ 
tain  an  IMAX  by  UMAX  grid.  The  material  in  the  rows  (or 
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columns)  to  be  added  will  have  the  same  properties  as  the  last 
row  (or  column)  of  the  grid.  Therefore  a  package  that  extends 
to  the  edge  of  the  grid  in  the  x-direction  before  rezone  is 
assumed  to  be  infinite  in  that  direction.  Likewise  for  the 
y-direction.  The  momentum  and  total  energy  of  the  original 
grid  are  conserved  when  the  cells  are  combined  by  converting 
some  kinetic  energy  into  internal  in  a  process  loosely  called 
''thermalizing . "  The  DX,  X,  TAU  and  TX  arrays  are  redefined 
when  the  grid  is  rezoned  in  the  x-direction.  For  the  y-direction 
rezone,  DY,  Y  and  TY  arrays  are  redefined. 

All  the  cell  quantity  arrays  (AMX,  AIX,  U,  V,  P,  MFLAG, 
STRSZZ,  STRSRR,  STRSRZ)  are  redefined  when  rezoning  in  one 
or  both  directions.  This  REZONE  routine  is  very  limited  and 
should  be  used  only  to  add  more  void  or  more  material  to  a 
semi-inf inite  package.  All  finite  boundaries  of  a  package 
must  be  within  the  original  grid. 

REZONE  is  called  from  EDIT  when  certain  input  parameters 
are  set  (NUMREZ,  REZFCT,  SS4) .  It  can  be  "forced"  on  the  first 
cycle  of  a  restart  run  by  setting  SS4=1.  Or  it  can  be  "triggered" 
by  PH2  which  senses  when  a  finite  signal  reaches  the  edge'of 
the  grid. 

7.1.23  SETUP 

This  routine  generates  problems  involving  one  spherical 

[21 

package  and/or  several  cylindrical  packages.  (The  "CLAM" '• 
generator  has  been  modified  to  set  up  HELP  problems  that  in¬ 
volve  more  complex  geometries.)  SETUP  also  defines  the  X, 

Y,  DX,  DY  and  TAU  arrays  for  variable  as  well  as  constant 
zoned  grids.  The  properties  of  each  cell  in  a  material  package 
are  determined  by  the  input  cards  read  by  SETUP.  The  coordinates 
of  the  material  tracers  are  defined  as  each  package  is  generated. 
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i 

Problems  which  generate  a  void  package  require  additional 
input  cards  to  define  the  "free  surface"  tracers.  However, 
when  the  void  surrounds  only  a  spherical  package,  these 
cards  are  unnecessary  as  the  free  surface  tracers  are  de¬ 
fined  when  the  sphere  is  the  first  package  generated.  Both 
SETUP  and  the  "CLAM"  generator  move  the  tracers  very  slightly 
off  of  the  grid  lines  to  simplify  the  logic  of  INFACE.  The, 
input  cards  read  by  SETUP  are  described  in  Section  7.2.3. 

7.1.24  STRNG 

The  yield  strength  of  the  material  in|a  cell  is  computed 
by  STRNG.  (The  strength  constants  for  each  material  are  de¬ 
fined  by  input  cards  when  the  problem  is  generated.)  The 
strength  of  material  in  a  mixed  cell  is  a  volume  weighted 
average  of  the  strength  of  the  materials  in  the  cell.  STRNG 
is  called  from  PH3.  , 
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7 . 2  GENERATING  PROBLEMS 

7.2.1  General  Capabilities  and  Limitations 


The  HELP  problem  generator  can  set  up  one  sphere  and 
*  *  * 

any  number  of  cylindrical  packages  or  subpackages  with  the 
following  restrictions. 


1.  All  cells  within  a  package  will  have  the 
same  properties. 

2.  The  right  and  top  grid  boundaries  cannot 
serve  as  package  boundaries.  If  a  package 
extends  to  the  right  edge  of  the  grid  the 
package  is  assumed  to  be  infinite  in  the 
x-direction.  Likewise  for  the  top  edge  of 
the  grid. 

3.  The  grid  .should  have  at  least  three  columns. 
(Special  conditions  for  1-D  problems  have 
not  been  included  in  the  code.) 

4.  The  right  and  top  grid  boundaries  are 

.  transmittive  although  the  bottom  boundary 
can  be  either  reflective  or  transmittive. 

5.  Void  opening  and  closing  mechanisms  have 
not  been  generalized  and  included  in  the 
current  version  of  HELP. 


6.  The  code  defines  only  one  y  "the  ideal 

gas  equation  of  state.  Problems  with  more 
than  one  ideal  gas  require  code  changes. 

It  might  be  noted  here,  also  that  the 


Limited  only  by  the  dimensions  of  the  mixed  cell  arrays. 

In  plane  coordinates  circular  and  rectangular  packages 
are  generated. 
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essentially  infinite  density  obtained  by 
iterating  for  pressure  equilibration  in 
a  cell  containing  cold  (E  =  0)  and  hot 
gas  may  lead  to  spurious  mass  transports. 
More  gradual  energy  changes  are  necessary 
for  gases.  ; 

7.  The  center  of  a  sphere  package  must  be 
on  a  horizontal  grid  line  and  on  the 
axis  of  symmetry.  Furthermore,  jif  the 
sphere  boundary  is  a  free  surface,  the 
sphere  must  be  the  first  package 
generated. 

8.  The  DATA  statements  in  BLOCK  and  EQST  must 
be  changed  if  the  user  wishes  to  redefine 
the  normal  density  or  the  equation  of  state 
constants  of  a  material. 

t 

The  following  capabilities  of  HELP’ should  also 
be  noted  here:  ' 

I 

1.  The  number,  of  material  packages  I  generated 
is  limited  only  by  the  dimensions  of  the 
variables,  which  can  be  increased  up  to 
the  limit  of  the  user's  computer  core 
storage. 

2.  A  cell  can  contain  any  number  of  materials 

and  material  interfaces.  ; 

3.  Two  packages  can  actually  contain  the 
same  material  (use  the  same  equation  of 
state)  but  have  different  initial  velo¬ 
cities,  energies  or  densities.  ■ 

i  -  . 
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4.  A  package  can  be  divided  into  any  number 

o£  spatially  disconnected  subpackages.  (See 
definition  of  MAT  in  the  following  discussion 
of  input  variables.) 

5.  When  strength  effects  are  unimportant,  the 
PH3  subroutine  can  be  bypassed  by  setting 
to  -1.  the  input  variable  CYCPH3.  In  this 
case  the  code  is  "hydrodynamic"  in  that  it 
solves  problems  in  compressible  fluid  dynamics. 

7.2.2  Format  of  Input  Cards 

A  HELP  problem  is  generated  by  three  kinds  of  input 

data; 

1.  The  identification  header  card  read  by  INPUT 
that  contains  alphanumeric  symbols  in  columns 
2  through  72 . 

2.  The  input  data  read  by  subroutine  CARDS  that 
defines  variables  in  the  Z-block  (the  first 
150  words  of  blank  common) . 

3.  The  input  data  read  by  subroutine  SETUP 
that  defines  the  DX,  DY  arrays  when  either 
or  both  are  variable,  the  boundaries  and 
properties  of  the  material  packages,  and 
the  coordinates  of  certain  free  surface 
tracer  particles. 

The  input  cards  read  by  CARDS  use  the  following  format: 

(II)  (15)  (II)  (E9./)  (E9.4) 

Col.  1  Cols.  2-6  Col.  7  Cols.  8-16  Cols.  17-35  etc 
L  M  N  Z(M)  Z(M+1) 
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The  value  o£  L  for  each  input  variable  is  listed  in  the  dis¬ 
cussion  that  follows.  When  L  =  0  the  variable  defined  is 
typed  real.  When  L  =  2,  the  variable  defined  is  typed  integer 
When  L  =  1  the  variable  defined  is  typed  real  and  the  card 
is  the  last  one  to  be  read  until  CARDS  is  called  again.  The 
value  of  M  is  the  location  in  blank  common  of  the  variable 
being  defined  in  Cols.  8-16;  this  location  number  is  also 
listed  below  for  each  input  variable.  The  value  of  N  indi¬ 
cates  how  many  consecutively  stored  variables  of  the  same  type 
are  defined  by  a  card  (usually  N  =  1) .  Z (M)  must  be  punched 
with  an  E  format  even  when  it  is  defining  an  integer  variable 

The  format  of  the  cards  read  by  SETUP  varies- and  is 
described  in  the  following  discussion  of  thei  input  variables. 

7.2.3  Description  of  Input  Variables 

The  following  list  of  input  variables;  is  ordered  as  the 
HELP  setup  input  deck  should  be  ordered.  The  first  group  of 
cards  are  read  by  CARDS,  the  second  group  by  SETUP  and  the  last 
dummy  card  is  read  by  CARDS.  Many  cards  in  the  first  group 
can  be  omitted  when  the  user  is  inputting  a  zero  value. 
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Input  Read  by  CARDS 


Variable 

Name 

Column 

One  Flag 

Location  in 
Blank  Common 

Definition 

PKCl) 

1 

151 

Problem  number.  Any  number  from  .0001  to 
99.9999.  Should  be  identical  to  PROB.  This 
is  the  second  input  card  of  a  restart  as  well 
as  a  setup  deck. 

PROB 

1 

Problem  number. 

NFRELP 

2 

S 

Gives  frequency  of  "long"  EDIT  prints  which 
gives  velocities,  energy,  compression  and 
stresses  for  all  cells  in  active  grid.  A 
"short"' EDIT  gives  this  information  only  for 
the  axis  column.  When  NFRELP=2,  every  second 
EDIT  print  is  "long". 

NDUMP7 

2 

6 

Gives  frequency  of  restart  tape  dumps. 

Program  will  dump  only  when  EDIT  prints. 

If  NDUMP7=5,  a  restart  dump  will  be  made 
every  Sth  time  EDIT  prints. 

ICSTOP 

2 

7  ■ 

Gives  cycle  at  which  execution  will  stop  if 
the  user  wants  to  specify  a  cycle  rather  than 
a  time  on  which  to  stop.  If  stopping  on  time 
omit  this  card. 

NUMREZ 

2 

12 

The  total  number  of  times  the  grid  will  be 
rezoned.  If  not  rezoning  the  grid  omit 
this  card. 

KUNITR 

2 

The  tape  unit  INPUT  will  read  from.  KUNITR 
will  be  set  equal  to  7  unless  defined  by  this 
input  card.  If  restarting  from  a  HLPCLM  tape, 
KUNITR  must  be  7. 

IPR' 

■  2  .  . 

IS 

Maximum  number  of  iterations  CDT  will  perform 
when  attempting  to  equilibrate  the  pressures 
of  materials  in  a  mixed  cell.  If  the  pressures 
are  not  within  an  epsilon  (PRCNT)  of  the 
average  pressure  after  IPR  iterations,  an  error 
exit  occurs.  (Try  3S.) 

PRCNT 

16 

Convergence  minimum  for  the  pressure  equili¬ 
bration  of  materials  in  a  mixed  cell.  All  pre 
pressures,  P. ,  must  satisfy  the  following: 
l(Pav  -  Pi)/Pavl  <  PRCNT.  (Try  .001.) 

KUNITW 

2 

17 

The  tape  unit  EDIT  and  SETUP  will  write  on. 
KUNITW  will  be  set  equal  to  7  unless  defined 
by  this  input  card.  If  restarting  from  a 

HLPCLM  tape  this  tape  unit  should  not  be  7 
unless  )the  user  wants  SETUP  to  write  over 
the  HLPCLM  dump. 

IGM 

2_ 

21 

When  IGM  =  0.  code  uses  cylindrical 
coordinates.  (Vhen  IGM  =  1.,  code  uses 
plane  coordinates. 

DMIN 

24 

The  maximum  relative  error  in  the  energy  sum 
that  the  user  wants  to  tolerate.  The  error 
is  computed  as  follows; 

('  KMAX  . 

-  EThJ/ETH 
K»1  / 

where  Ew  is  total  energy  of  cell  K  and 
ETII  is  theoretical  energy  in  >the  grid- -com¬ 
puted  in  SETUP  and  updated  in  PUl,  PH2,  PH3. 
Usually  DMIN  -  10'*. 
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Variable 

Name 

Column 

One  Flag 

Location  in 
Blank  Common 

Definition 

CVIS 

27 

A  flag  that  describes  the  boundary  condition  at 
the  bottom  of  the  grid.  If  CVIS  ■  0.,  the 
boundary  will  be  reflective.  If  CVIS  ■=  -1., 
the  boundary  will  be:  transmittive . 

IMAX 

2 

33 

The  number  of  columns  in  the  grid.  If 
planning  to. rezone  the  grid  once  IMAX  should 
be  an  even  number;  IMAX,  for  this  version 
of  the  program,  should  be  at  least  3. 

JMAX 

2 

3S 

1 

The  number  of  rows  in  the  grid.  (Follow  same 
rules  as  for  IMAX  if  grid  is  to  be  rezoned. ) 

MAPS 

2 

42 

When  MAPS  >0.,  part  of  the  EDIT  print  is  a 
set  of  symbolic  maps  of  the  compression, 
pressures,  velocities,  and  specified  internal 
energies  of  cells  in  the  active  grid. 

NUMSCA 

2 

43 

The  number  of  times  jthe  cycle  or  time  interval 
between  EDIT  prints  is  increased.  (See  PRLIM) 

PRLIM 

44 

The  time  or  cycle  at  which  the  EDIT  print  fre- 

quency  is  changed.  Example:  you  wish  to  print 
every  10"*  sec  until  T  =  10'^  sec;  and  every 
10*’  sec  until  T  =  10*‘  sec,  and  every  10"*sec 
thereafter.  I 

Set:  PRDELT  =>  10' * 

IPCYCL  “0. 

PRLIM  -10*’ 

PRFACT  -10. 

NUMSCA  “2. 

NOTE:  When  PRDELT  is  increased  by  a  factor, 
PRLIM  is  also  increased  by  the  same  factor  ■ 
for  the  next  rescaling.  If  you  want  a  con¬ 
stant  print  interval  omit  NIIMRCA,  PRLIM, 
PRFACT. 


PRDELT 

45 

The  time  (sec)  between  EDIT  prints.  If 
printing  on  cycle  intervals  omit  this  card. 

PRFACT 

46 

The  factor  by  which  the  print  interval 
(PRDELT  or  IPCYCL)  and  the  print  limit 
(PRLIM)  are  increased.  (See  PRLIM) 

11 

2 

47 

The  number  of  columns  in  the  grid  that  have 
non-zero  velocities  of  energies  plus  2.  11 

and  12  define  the  "active"  grid.  (Most  cal¬ 
culations  are  done  only  inside  the  active 
grid.) 

12 

2 

48 

The  number  of  rows  in  the  grid  that  have 
non-zero  velocities:  or  energy  plus  2^. 

IPCYCL 

2 

49 

The  number  of  cycles  between  EDIT  prints. 

If  printing  on  time  intervals  omit  this 
card. 

IVARDY 

2 

54 

When  the  grid  is  variably  zoned  in  the  y- 
direction  (DY's  not  constant) ,  IVARDY  =  1. 
When  the  grid  has  constant  zoning  in  the 
y-direction,  this  card  is  omitted. 

(IVARDY  -  0.) 

98 


Variable 

Name 


Location  In 
Blank  Common 


Definition 


3SR-350 


Column 
One  Flag 


N6 

2 

S6 

This  parameter  causes  the  program  to  set 
negative  pressures  to  zero  in  the  rows  of  the 
grid  below  J  »  N6.  Omit  this  card  if  you 
want  to  allow  negative  pressures  everywhere 
in  the  grid.  If  you  want  to  set  negative 
pressures  to  zero  everywhere,  set  N6  •  10*. 

GAMMA 

-- 

62 

(a+1.)  in  gamma-law  equation 
of  state.  Used  only  for  an 
ideal  gas. 

NSIDES 

2 

67 

The  number  of  line  segments  that  describe  the 
void  package  not  counting  the  grid  boundary 
lines.  Omit  this  card  if  there  is  no  void  or 
if  the  only  free  surface  is  a  spherical 
surface. 

NMAT 

2 

68 

The  number  of  material  packages,  excluding  the 
void  package.  The  maximum  number  of  packages 
is  determined  by  the  dimensions  of  the  mixed 
cell  variables  in  MCXELL  COhlUON  and  in  BLANX 
COMMON.  The  minimum  number  is  two. 

CYCMX 

69 

The  number  of  passes  through  subroutine  INFACE 
in  one  time  step  to  minimize  transport  noise 
near  interfaces.  Usually  CYCMX  =  2,  but  when 
calculation  costs  permit,  CYCMX  =  8  is  a 
good  maximum. 

CYCPH3 

70 

The  number  of  passes  through  subroutine  PH3 
(strength  phase)  in  one  time  step.  Usually- 
CYCP113  =  1.,  but  larger  values  can  help  to 
dampen  oscillations.  When  the  calculation  is 
pure  hydro  (no  strength  effects)  set  CYCPH3  = 
-1. 

REZFCT 

71 

A  rezone  flag.  When  the  grid  is  to  be  rezoned 
at  some  point  in  the  calculation,  REZFCT  =  1, 
Omit  this  card  only  if  the  grid  will  never  be 
rezoned. 

NTRACR 

2 

72 

The  number  of  tracer  particles  per  cell  edge 
to  be  used  to  circumscribe  the  material 
packages.  Include  this  card  even  when 

the  problem  is  generated  by  HLPCLM,  since 
this  parameter  is  used  in  adding  tracers 
(ADDTCR)  as  well  as  in  setting  up 
material  packages  (SETUP) .  (Try  2  to 

S.) 

NMXCLS 

2 

73 

The  maximum  number  of  mixed  cells  to  exist  in 
the  grid  at  any  one  time  during  the  calcula¬ 
tion.  Consult  the  dimensions  of  the  mixed 
cell  arrays  in  MXCELL  COMMON  and  BLANK  COMMON. 

NTPMX 

2 

78 

The  maximum  number  of  tracer  particles  per 
material  package.  This  maximum  should  depend 
upon  the  dimensions  of  the  TX  and  TY  arrays. 

GLUED 

79 

When  GLUED  =  1.  all  mixed  cells  are  "glued" 
after  PHI  updates  velocities  and  energies 
due  to  pressure  gradients.  This  card  is 
usually  omitted  since  the  program  will  auto¬ 
matically  glue  underdense  cells  that  are 
"kicked”  by  PHI.  This  option  has  been  used 
when  all  the  mixed  cells  were  between  an 
expanding  hot  gas  and  a  cold  dense  material. 

NTCC 

2 

81 

When  NTCC  •  1.,  passive  cell  centered  tracer 

particles  are  initially  placed  in  the  center 
of  every  fourth  non-empty  cell  and  subsequent¬ 
ly  moved  across  the  grid  with  the  material. 
These  are  used  only  for  producing  particle 
plots.  When  this  card  is  omitted  (NTCC  -  0), 
cell  centered  tracers  are  not  computed. 
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Variable 

Name 

IVARDX 

EMIN 

PMIN 

INTER 

NODUMP 

ROEPS 

FINAL 

NADO 

MINX 

MAXX 

MINY 

MAXY 

lEXTX 

JEXTY 

SS4 

DXF 

DYF 

STAB 


Column  Location  in 

One  Flag  Blank  Common  Definition 


2  83 

85 

86 

2  87 


2  96 

110 

113 

2  118 


2 

2 

2 

2 


119 

120 
121 
122 


2  123 


2  124 


130 


136 


137 


When  the  grid  is  variably  zoned  in  the  x- 
direction  (DX's  not  constant),  IVARDX  -  1. 

When  the  grid  has  constant  zoning  in  the 
x-direction,  this  card  is  omitted  (IVARDX  =  0.). 

The  minimum  specific  internal  energy  of  an 
ideal  gas  to  be  used  in  the  pressure  iteration. 
Usually  EMIN  =10’. 

A  pressure  cutoff  for:  the  zero  cycle  of  the 
calculation.  Thereafter  PMIN  is  computed  in 
CDT.  Usually  PMIN  =  10*. 

A  special  editing  flag  for  debugging.  When 
INTER  =1..  EDIT  prints  after  PHI  and  PH3 
as  well  as  after  CDT  on  print  cycles.  When 
INTER  =  99.,  PH3  prints  energy  totals  after 
updating  each  cell;  PH2  does  the  same  when 
INTER  »  7.  The  later  options  are  in  addition 
to  the  extra  EDIT  prints  and  involve  many 
pages  of  output. 

When  NODUMP  =  1.,  EDIT  will  not  write  any 
tape  dumps.  (This  flag  overrides  the  NDUMP7 
option,  but  does  not  prevent  SETUP  from 
writing  the  cycle  0  dump.) 

A  round-off  epsilon  used  to  define  cutoffs 
in  the  calculation.  Try  ROEPS  =  10'*. 

The  final  value  of  the  stability  fraction 
used  in  determining  the  time  step.  Usually 
FINAL  -  .4.  (See  STAB) 

A  flag  for  adding  tracer  particles.  If 
NADO  »  3,  EDIT  calls  ADDTCR  on  the  first  cycle 
of  - the  run,  ADDTCR  add,s  3  tracers  between  all 
existing  pairs  of  tracers  in  a  specified 
region  (except  where  there  will  he  more  than 
NTRACR  tracers  in  a  cell.)  NADD  is  reset  to 
zero  in  ADDTCR. 


These  parameters  specify  the  left  and  right 
columns  and  the  bottom  and  top  rows,  respec¬ 
tively,  of  the  region  in  which  ADDTCR  will 
add  tracer  particles.  This  card  must  be 
included  when  NADD  ^0. 

A  rezone  flag.  The  grid  will  be  rezoned  in 
the  x-direction  only  iihen  lEXTX  =  1.  Omit  this 
card  if  the  grid  is  not  to  be  rezoned  or  if  it 
is  to  be  rezoned  in  the  y-direction  only. 

A  rezone  flag.  The  grid  will  be  rezoned  in  the 
y-direction  only  when  JEXTY  =  1.  Omit  this 
card  if  the  grid  is  not  to  be  rezoned  or  if  it 
is  to  be  rezoned  in  the  x-direction  only.  . 

A  rezone  flag.  The  grid  is  rezoned  on  the 
first  cycle  of  the  run  when  SS4  =  1.  (This 
overrides  REZFCT  but  not  NUMREZ  which  must 
be  non-zero  for  REZONE  to  be  called.) 

The  x-dimcnsion  of  cells  (cm)  when  the  grid 
has  constant  zoning  in  the  x-direction.  Omit 
this  card  if  the  x-dimension  of  cells  varies. 

The 'y-dimension  of  cells  (cm) 
when  the  grid  had  constant 
zoning  in  the  y-direction .  Omit 
this  card  if  the  y-dimension  of 
cells  varies. 
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"stability  fraction"  for  the 
calculation  of  DT.  If  FINAL=0.,, 

STAB  is  constant.  Otherwise  its 
value  increases  from  fhe  input 
100  to  FINAL  in  a  geometric  progression. 
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Variable 

Name 

Column 

One  Flag 

Location  in 
Blank  Common 

Definition 

DTMIN 

144 

The  minimum  value  for  the  time  step.  The 
program  gives  an  error  exit  if  CDT  calculates 
a  At  <  DTMIN.  (Try  DTMIN  =  10'“) 

JCEKIR 

2 

147 

The  J-line  that  passes  through  the  center  of 

the  sphere  package.  (Used  in  GLUE  and  PH2 
to  insure  symmetric  results  in  spherical  cal- 


culations.)  Omit  this  card  if  the  problem 
does  not  involve  a  spherical  package. 

RADIUS 

148 

Radius  (cm)  of  spherical  package.  Omit  this 
card  if  the  problem  does  not  involve  a 
spherical  package. 

BBAR 

149 

A  constant  used  in  calculation  of  local 
sound  speed  of  all  materials  except  ideal 
gases.  C  =  C  +  BBAR  •  (/F) .  (Usually 

BBAR  =0.5)  “ 

TSTOP 

1 

50 

The  time  (sec)  at  which  the  calculation  will 
stop.  This  card  (with  a  "1"  in  column  one) 
should  not  be  omitted  even  when  stopping  on 
a  specified  cycle,  and  should  always  be  the 
last  card  of  a  restart  deck. 

(Insert  all  cards 

described  on 

the  following  pages) 

EMOB 

1 

150 

Dummy  end  card  of  SETUP  deck.  Omit  from 
restart  deck.  EMOB  =  0. 
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Input  Read  by  SETUP 

1.  Defining  Cell  Dimensions  When  Variable.  User 
can  specify  variable  DY  and  /or  DX  values.  In  the  case  that 
both  are  variable,  the  ^  array  is  defined  first.  Both 
arrays  are  read  using  the  following  format:  ; 


Variable 

Name  Columns  Definition 


NT(I),  1=1,4 

1-4 

5-8 

NT(I) 

is  the  number  of  zones  that  have  dimension 

.  TEMP(I).  When  all  DY's  are  defined,  set  next 

9-12 

999.  Likewise  for  DX's, 

13-16 

1  NT  = 

1 

TEMP(I),  1-1,4 

17-30  1 

NOTE : 

When  DY's  are  variable,  exactly  JMAX  DY's 

31-40 

must 

be  defined.  Likewise,  exactly  IMAX  DX's 

41-50 

must 

be  defined. 

51-60 

2.  Defining  Each  Material  Package  (4  cards  for 


each  package  or  subpackage) . 


Variable  Card 

Name  Number  Columns  Definition 


IGEOM 

1 

1 

MAT(N) 

2 

1-6 

PACRT 

2  7-16 

PACTP 

17-26 

PACLF 

27-36 

PACBT 

37-46 

Defines  package  geometry.  IGEOM  =  1  for  a  rectangle. 
IGEOM  =  2  for  a  sphere  (II). 

Material  code  number  of  this  package  (e.g.,  if 
|MAT(1) 1  =  3  then  package  one  is  iron.  See  list 
in  EQST) .  If  material  package  N  is  divided 
into  several  spatially  disconnected  subpackages, 

MAT(N)  is  negative  for  all  but  the  last  subpackage 
(16). 

(Omit  for  spherical  package.)!  Location  (cm)  of 
right,  top,  left,  bottom  boundaries,  respecitvely , 
of  rectangular  package.  (4E10.4) 

Set  PACRT  =  -1.  if  package  is  infinite  in  x-direction. 
Set  PACTP  >=  -1.  if  package  is  infinite  in  y  direction. 
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Variable  Card 

Name  Number  Columns 


Definition 


UUR(N) 

3 

1-10 

Radial  velocity  of  all  material  in  the  package. 

(E10.4) 

VVA(N) 

11-20 

Axial  velocity  of  all  material  in  the  package. 

(E10.4) 

SSIEN(N) 

21-30 

Specific  internal  energy  of  all  material  in  the 
package  .(E10.4) 

RHOIN(N) 

31-40 

Initial  density  of  al,l  material  in  the  package. 

(E10.4) 

CZERO(N) 

4 

1-10 

Y  . 

0 

Yield  strength  of  all  material  in  package 

N  - 

STKI(N) 

11-20 

Y 

1 

(y  +  Y  u  +  Y  pM  (l  -  E/E  ) 

STK2(N) 

21-30 

Y 

i 

■  where 

STEZ(N) 

31-40 

E 

0 

p  =  p/p  -  I. 

0 

■ 

E  “  specific  internal  energy  .C4E10 . 4) 

RMUCN) 

41-SO 

Rigidity  modulus  of  all  material  in  package  N. 
(E10.4) 

AMDMCN) 

51-60 

Failure  criterion  based  on  compression  for  all 

material  in  package  N.  If  compression  of  material 
in  a  cell  is  less  than  AMDM,  the  cell  pressure 
cannot  be  negative.  (Usually  AMDM  =  0.9  to  0.99). 
CE10.4) 


3.  Defining  Tracer  Particles  Along  Each  Line  Segment 
of  the  Void  Package  Boundary.  (One  card  for  each  line  segment 
Omit  these  cards  when  NSIDES  =  0.) 


Variable 

Name  Columns  Definition 


STRTX 

1-10 

x-coordinate 
STRTX  =  -1. 

(cm) 

of 

starting 

point . 

If 

point  is  at  infinity 

STRTY 

11-20 

y-coordinate 
STRTY  =  -1. 

(cm) 

of 

starting 

point. 

If 

point  is  at  infinity 

ENDX 

21-30 

x-coordinate 
ENDX  =  -1. 

(cm) 

of 

end  point 

.  If 

point 

is  at  infinity, 

ENDY 

31-40 

y-coordinate 
ENDY  =  -1. 

(cm) 

of 

end  point 

.  If 

point 

is  at  infinity. 

(1)  SETUP  will  read . "NSIDES"  cards,  one  for  each  line  segment  of  the  free  surface. 

(2)  The  input  of  these  tracers  should  be  ordered  such  that  when  one  travels  between 
consecutive  points,  the  void  is  on  the  left. 

(3)  If  the  void  boundary  is  spherical  the  tracers  for  the  void  package  are  automatically 

defined  the  spherical  package  is  generated. 
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7.2.4  Sample  Input  Deck  for  Setting  up  a  HELP  Problem 


The  input  for  setting  up  an  impact  of  an  iron  cylinder 
onto  a  one-centimeter  aluminum  plate  is  listed  below.  This 
impact  problem  has  been  run  on  the  version  of;  the  code  listed 


in  this  report  and  is  a  suggested  test  problem  for  the  new 
user  of  the  HELP  code.  ' 


1 

151 

126.5 

1 

126.5 

5 

12. 

2 

6 

12. 

2 

15 

140. 

16 

t.oni 

?ii 

11.0  -0  3 

27 

1-1. 

2 

33 

131. 

2 

35 

162. 

2 

42 

11. 

45 

15.0  -07 

2 

47 

7. 

2 

44 

130. 

2 

67 

IS. 

2 

63 

12. 

65 

12. 

70 

11. 

2 

72 

12. 

2 

73 

1250. 

2 

78 

1400. 

2 

81 

11. 

2 

83 

1  . 

46 

1.0  *06 

lie 

1.0  -05 

1131 

.4 

137 

.05 

1351 

i.n  -03 

1441 

1.0  -11 

1451 

.5 

1 

5t]|] 

2.0  -06 

TRON  ALUMINUM  IMPACT 


> 


20 

1  1 

1 

7.5 

-024.6 

-029.9 

-021  .14 

-01 

t 

1  1 

1 

1.31 

-011.51 

-011.73 

-Cl  1.99 

-01 

1 

599 

1  1 

1 

2.25 

-012.64 

-013  .04 

rOll  .50 

-01 

1 

3 

• 

37  5 

1  .40 

0. 

.  25 

! 

0. 

5.0 

♦  040. 

7.8 

6  .0 

1 

♦  0911. 

0. 

3.0 

♦  101.93 

412.97 

4 

3 

.25 

2.4 

C. 

1.  4 

0. 

0. 

0. 

2.79 

, 

3.0 

♦090. 

0. 

7.0 

♦  092.74 

Vll.985 

0. 

2.4 

3.25 

2.4 

3.25 

2.4 

3.25 

1.4 

! 

3.25 

1.4 

.375 

1  .4 

.375 

‘  1.4 

.375 

.25 

< 

.375 

.25 

0. 

.2  5 

1 

Ij  15010. 
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The  lines  that  follow  are  part  of  the  printed  output 
from  the  sample  problem  described  above.  This  calculation 
was  made  on  a  UNIVAC  1108. 
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PPCeLE” 

2A.S)Cn"‘ 

T  1  '-F 

.oonc2r :-p7 

f.YfiF  rrT.EN-THnOR.  ,,A}<  •RrL,tRhnN-CVCL» 

VE,  1  /♦  1  0  b.'TO  1  7bNA-n7  N 

It  j  L  t  ’  6  - » 

-b,9*;H;i26ti*02 

•h^  t  L  a  •.  !  U  ri  as  t  I  C  .0»<k 

l,3H|7b2*->»UH 

p^ckagf  no, 

\ 

7 

i  E 

3f  IbCSaip^Ul 
b.rif^6'M7*Ub 

kf 

1  ;v7?^9^j3 

tot, IN.  (SUMJ 
\  .•<9(,w*.Y?2»10 

1  •  1  NI,  1  1  90*U9 

•'ASS 

3*9b7H2JC1*C0 

9.2S90bVV*C* 

nv 

j,3t,ti9/2N*Ub 
2*07SA7 1 3*0H 

;u  J  Pj'S  n  1  Vt ) 
J*3Sb772-4*UL> 

2  *  1  N209OS  +  0N 

f.U 

1  ,  Llljbi  (j2  1  *0** 

1  ,365236e+0H 

PuASI  lC-*C-iA 
S,./PA':V32».jf 
7  t JSSnnS 1 »a7 

TOTALS 

8, 36''2257»0G 

I  .S2]2r,ip.i:, 

1 

9.».S‘43H22»rl 

3.b/A<  '.439S*0S 

3.S7jlfll‘*tCS 

2 , 37  1  M392*0‘« 

1*3090'436«08 

BCllf.'OAPV 

bcttoh 

right 

TOP 

SE  VAPnRATt  :}^ 

HASS  out 
ENERGY  OUT 
HU  OUT 

HV  OUT 

O.OOCLODC 

Q ,  OOCCOOO 
040000000 
0,0000000 

c.coocrnc 

C*CCOOOOC 

c • ocoorco 
c^oocDcro 

c  .occccor. 
0*0000000 
OtOOCOCOC 
0*0000000 

r .  conunr.o 

0  .rocL-oL'n 
O.DOOOOOO 
o.GOConoo 

ftORK  OO^^C 

o«oocoooo 

C*000GC00 

0*0000000 

T*Pt  7  DUMP  ON  CTCLT  2£. 


CELL-COCKoINaTCS  Of  TPACtRS  FOK  EACH  haTERIAL  PaCKa(iE 


PACKAGE  1 


N 

X 

T 

N 

X 

1 

•  UD 

S  *  90 

7 

.SO 

6 

2.60 

6.9C 

7 

3-00 

1  1 

S.OC 

S  .90 

12 

5.00 

16 

■  5.00 

'SiMO 

17 

-5.00 

21 

S.OC 

10*90 

22 

5.00 

26 

5 .00 

13. MO 

27 

5.00 

3  1 

S.OC 

1S.90 

32 

5.00 

36 

b  *00 

1  H.MC 

37 

5.00 

H  1 

.  s.oo 

20.90 

M2 

5.00 

H6 

S.DO 

23, 3f 

M7 

b.OO 

SI 

S.OM 

2S.R3 

S2 

5.0M 

S6 

S*i9 

2B.  10 

57 

5.20 

61 

3*03 

28.60 

62 

2.53 

66 

•  SO 

28.6  t 

67 

•  oo 

Y 

N 

X 

Y 

U 

5.9C 

3 

1  .00 

5 .911 

H 

5,90 

6 

3.50 

S.90 

V 

6. MO 

13 

5.00 

6.90 

M 

0.90 

18 

5.00 

-  9  .  MO 

19 

1  1  •  MC 

23 

S.OO 

1  1.90 

2H 

13.90 

28 

5.00 

IM.MO  ' 

2V 

16.M0 

33 

5.00 

16.90 

3M 

IB. 90 

38 

S.OO 

19. MO 

39 

21  *90 

M3 

s.oo 

21.69 

MM 

23.88 

M8 

s.oo 

2M.37 

M9 

26.27 

S3 

s.os 

26.7  1 

SM 

28. Ml 

58 

M.7I 

20. M6 

59 

20.62 

63 

2.02 

28.61 

6M 

28.6t 

66 

•  00 

S.90 

X 

i 

N 

X 

Y 

1  .50 

b.90 

5 

2*00 

S»9U 

M.DO 

5*90 

ID 

M.5D 

S*9D 

5.00 

7  »MU 

lb 

5*00 

7 .90 

5.00 

9.90 

20 

5.00 

- lO.MO 

5.00 

1  2  .  MO 

25 

5*00 

12*90 

5.00 

IM.90 

30 

5*00 

I5*MU 

5.00 

17.H0 

3S 

StCQ 

17*90 

b  .UD 

19.90 

MO 

5*00 

20*90 

5*00 

22*39 

Mb 

5.00 

22*«9 

5.01 

2M*86 

50 

5*02 

25*35 

5. 1 1 

27*2U 

55 

5*17 

27*70 

M.I3 

28*52 

60 

3*53 

28*58 

1.50 

28*61 

65 

1*00 

28*61 
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7.3  RESTARTING  PROBLEMS 


During  the  calculation,  HELP  periodically  writes  on  a 
tape  or  drum  file  the  problem  parameters  and  the  current 
state  of  the  material  in  each  cell.  By  reading  this  file 
and  a  few  input  cards,  the  code  can  "restart"  and  continue  a 
calculation  from  an  intermediate  point. 

The  three  cards  that  must  be  included  in  a  restart 
deck  are  as  follows: 

1.  the  heading  card 

2.  the  "restart"  card 

3.  the  last  card 

The  heading  card  can  contain  any  alphanumeric  symbols 
between  columns  2  and  72.  The  "restart"  card  defines  the 
four  variables  described  below. 


Variable 

Name 

Columns 

Definition 

PK(1) 

8-16 

Problem 

number. 

PK(2) 

17-25 

Restart 

cycle  number. 

PK(3) 

26-34 

Restart 

flag. 

=  -1.  prints  long  edit  of  restart  cycle. 

=  -2.  prints  short  edit  of  restart  cycle. 

=  -3.  restarts  from  HLPCLM  output  and 

prints  long  edit  of  restart  cycle. 

I 

PK(4)  35-43  Number  of  material  packages  in  the  problem 

(excluding  the  void  package) . 
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This  restart  card  is  read  by  CARDS  and  columns  1  through  7 
are  punched  as  follows  !• 

1234567  (Columns) 

1  1  5  1  4  ^ 

which  indicate  there  are  four  words  defined  by  the  card,  the 
first  one  stored  in  the  151st  location  of  blank  common. 

The  last  card  defines  TSTOP,  the  time,  to  stop  the  cal¬ 
culation.  This  must  be  the  last  card  of  the'  restart  deck 
even  when  the  problem  is  stopped  on  a  specified  cycle.  The 
TSTOP  card  is  punched  as  described  above  for’  a  setup  deck, 

e.g.  , 

1  2  3  4  5  6  7  8  9  10  12  13  14  15  16  (Columns) 

1  5011.0  -  06 

Many  problem  parameters  stored  in  blank  common  can  be 
redefined  when  restarting  at  an  intermediate  point  in  the  cal^ 
culation.  For  example,  the  number  of  times  PH3  is  subcycled, 
the  frequency  of  printing,  or  the  round-off  epsilon  can  be 
changed  by  including  in  the  restart  deck  cards  that  redefine 
these  variables.  These  extra  cards  must  be  added  between  the 
"restart”  card  and  the  last  card,  i.e.,  between  the  two  cards 
that  have  a  "1"  in  column  1.  All  cards  except  the  heading 
card  of  the  restart  deck  are  read  by  the  CARDS  subroutine  with 
the  format  described  in  Section  7.2.2. 

The  three  cards  listed  below  could  be  used  to  restart 
on  cycle  78  the  sample  problem  described  in  Section  7.2.4. 


Execution  of 

that 

problem 

would  continue  until  T  =  3.0  ysec. 

IRON 

ALUMINUM  IMPACT 

1 

151 

4j26.5 

78.  -2. 

1 

50 

13.0 

-  .06 
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7.4  DICTIONARY  OF  IMPORTANT  VARIABLES 


This  section  includes  a  description  o£  the  use,  units, 
dimension,  and  location  of  all  variables  in  named  or  blank 
commons  and  some  variables  local  to  the  subroutines.  The 
following  conventions  are  used  in  the  dictionary: 


(NAME) 

NAME 

B.C. 


=  ZCN} 


Variable  is  local  to  subroutine 
NAME. 

Variable  is  in  common  block  NAME. 

The  variable  is  in  Blank  Common 
or  equivalenced  to  a  variable  in 
Blank  Common. 

The  variable  is  equivalenced  to  a 
member  of  the  Z-array,  the  first 
array  in  Blank  Common.  These 
variables  are  usually  used  in 
setting  up  and  restarting  problems. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

AIX 

B.C. 

2500 

ergs/g 

Specific  internals  energy  in  a  cell. 

(IMAX  by  JMAX  array.) 

f  ZH 

ALE 

(MAP) 

41 

This  array  has  alphabetic  characters 
for  positive  values  in  the  density, 
velocity  and  energy  maps.  Defined 
in  a  DATA  statement. 

AMDM 

B.C. 

3 

INPUT  parameter.  .  A  material  i  with 
compression  >  AMDM.  is  considered 
•  solid.  i 

Usual  value:  0.95  to  0.99. 

Used  in  EQST  in  testing  whether  to 
allow  negative  pressures  (tensions) 
in  pure  cells.  Used  in  PH3  to  cal¬ 
culate  SOLID  =  AMDM  *  RHOZ,  which 
determines  whether  material  stresses 
are  nonzero. 

AMMP 

(PH2) 

1 

g 

Mass  transported  across  right 
boundary  of  a  cell.  (See  Appendix 

B),. 

AMMU 

(PH2) 

1 

cm-g/sec 

Radial  momentum  transported  across 
the  bottom  boundary  of  a  cell.  (See 

Appendix  B) . 

AMMV 

(PH2) 

1 

cm-g/sec 

Axial  momentum  transported  across 
the  bottom  boundary  of  a  cell  (See 

Appendix  B) .  . 

AMPY 

(PH2 

1 

g 

Mass  transported  across  top  boundary 
of  a  cell.  (See  Appendix  B) 

AMUR 

CPH2) 

1 

cm-g/sec 

Radial  momentum  transported  across 
right  boundary  of  a  cell.  (See 

Appendix  B)  J 

AMUT 

CPH2 

1 

cm-g/sec 

Radial  momentum  transported  across 
top  boundary  of  a  cell.  (See 

Appendix  B) 

AMVR 

(PH2) 

1 

cm-g/ sec 

Axial  momentum  transported  across 
right  boundary  of'  a  cell.  (See 

Appendix  B) 

AMVT 

CPH2) 

1 

cm-g/sec 

Axial  momentum  transported  across 
top  boundary  of  a  cell.  (See 

Appendix  B) 

AMX 

B.C. 

2500 

g 

Total  mass  in  a  cell.  (IMAX  by 

JMAX  array) 

BEAR 

-Z(149) 

1 

Used  in  CDT.  An  INPUT  parameter 
used  in  local  sound-speed  calculation 
for  all  materials  other  than  a  poly¬ 
tropic  gas.  (Local  sound-speed  for 
material  i  in  cell  K  is  approximated 
as  (Cgj^)  +  BEAR  •  /P  (Kj ,  where  the 
coefficient  B  is  obtained  by  de¬ 
termining  a  typical  slope  for  the 
isentropes  in  Refl  5  and  using  the 
relation  C  =  V/^dP/dV  to  evaluate 

5  at  a  particular  pointj 
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Variable  Number  of  ' 

Name  Location  Words  Units 


Definition 


BBOUND 

■Z(74) 

1 

ergs 

Calculated  in  PH3.  Printed  in  EDIT 
under  "Elastic  Plastic  Work."  Total 
work  done  by  the  elastic  and  plastic 
stresses. 

BOTH 

-2(39) 

1 

e 

Calculated  in  PH2.  Printed  in  EDIT. 
Total  mass  transported  across  bottom 
of  grid'. 

BOTMU 

»Z(64) 

1 

cra-g/sec 

Calculated  in  PH2.  Printed  in  EDIT. 
Total  radial-momentum  transported 
across  bottom  of  grid. 

BOTMV 

-2(40) 

1 

cm-g/sec 

Calculated  in  PH2.  Printed  in  EDIT. 
Total  axial -momentum  transported 
across  bottom  of  grid. 

CNAUT 

MXCELL 

30 

cm/sec 

Used  in  CDT.  Approximate  sound- 
speed  of  nineteen  materials  defined 
in  a  DATA  Statement  in  BLOCK. 

1  ESCAPA^  _ _ 

'^Oi  "  S  °  • 

CVIS 

=2(27) 

1 

INPUT  parameter.  Used  to  describe 
the  bottom  boundary-condition  in 

PHI,  PH2  and  PH3.  Bottom  boundary 
is  transmittive  when  CVIS  =  -1., 
reflective  when  CVIS  =  0. 

CYC 

B.C. 

1 

Defined  in  INFACE.  Used  in  FLUX  to 
define  flux  terms  of  cells  that  be¬ 
come  mixed  after  first  subcycle  of 
INFACE. 

CYCLE 

-2(2) 

1 

•• 

Used  in  INPUT,  SETUP,  CDT,  EDIT. 

Cycle  number  (an  integer  value  in 
floating  point  fcm) . 

CYCMX 

-Z(69) 

1 

INPUT  parameter.  Equals  number  of 
passes  through  INFACE  on  each  cycle. 
Suggested  minimum  of  2  and  maximum  of 
8.  Used  to  minimite  transport  noise 
near  interfaces. 

CYCPH3 

-2(70) 

1 

-- 

INPUT  parameter.  Number  of  times 
to  subcycle  PH3.  If  value  if  -1., 

PH3  is  omitted. 

CZERO 

B.C. 

3 

2 

dynes /cm 

Value  of  Yq  for  yield  strength 
calculation;  Defined  by  input  cards 
for  each  material  package  in  the  grid 
Used  in  STRNG.  (See  STRENG) 

DDX 

B.C. 

52 

cm 

An  array  equivalenced  to  the  DX  array 
--such  that  DDX(l)  =  DX(0). 

DDY 

B.C. 

102 

cm 

An  array  equivalenced  to  the  DY  array 
such  that  DDY(l)  =  DY(0). 

DELEB 

(PH2) 

1 

ergs 

Total  energy  associated  with  mass 
transported  across  bottom  boundary 
of  a  cell.  (See  Appendix  B) 

DELER 

(PH2) 

1 

ergs 

Total  energy  associated  with  mass 

transported  across  right  boundary 
of  a  cell.  (See  Appendix  B) 
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Variable 

Number  of 

Name 

Location 

Words 

Units 

Definition 

DELET 

(PH2) 

1 

ergs 

Total  energy  associated  with  mass 
transported  acro.ss  top  boundary  of 
cell.  (See  Appendix  B) 

DELI 

(PH2,  PH3) 

1 

ergs/g 

Change  of  specific  internal  energy 
of  a  cell.  I 

DELM 

(PH2) 

1 

g 

Total  change  in  mass  of  a  cell. 

DELU 

(PH2,  PH3) 

1 

cm/sec 

Change  of  radial  velocity  of  a 
cell. 

DELV 

(PH 2,  PH3) 

1 

cm/sec 

Change  of  axial  velocity  of  a 
cell. 

DMIN 

=  2(24) 

1 

■■ 

INPUT  parameter.  Allowable  relative 
error  in  energy  sum.  When  relative 
error  is  >  DMIN,  calculation  is  ter¬ 
minated.  The  energy  error  is  com¬ 

puted  in  EDIT.  !If  everything  is 
working  right  you  should  be  able  to 
use  10" ®  for  DMIN. 

DT 

=  2(3) 

1 

sec 

Time  step.  Calculated  in  CDT . 

DTFACT 

(PH3) 

1 

.  - 

Factor  used  in  calculating  a  variable 

time  step  when  subcycling  the  PH3 
calculations. 

DTMIN 

=  2(144) 

1 

sec 

INPUT  parameter;  Minimum  time  step; 
used  in  CDT.  After  STAB  =  FINAL, 
if  DT  <  DTMIN  execution  is  stopped. 

DTNA 

=  2(26) 

1 

sec 

DT  from  previous  time  cycle.  Used 
in  INPUT,  CDT,  and  EDIT. 

DTNOW 

(EDIT) 

1 

sec 

Used  for  saving  DT  when  calling  CDT 

to  recalculate  pressures  after  a 
REZONE. 

DTSTR 

Local 

1 

sec 

Used  in  PH3.  Time  step  for  a  PHI 
subcycle. 

DX 

B.C. 

SO 

cm 

The  radial-dimension  of  cells. 

Equivalenced  to  DDX  such  that  DDX(l) 

=  DX(0)  . 

DXF 

=2(136), 

1 

cm 

An  INPUT  parameter  used  to  calculate 
the  DX  array  if  the  radial  dimension 
of  the  cells  is  uniform. 

DY 

B.C. 

100 

cm 

The  axial -dimension  of  cells. 
Equivalenced  to  DDY  so  that  DnY(l) 

=  Dy(0). 

DYF 

=  2(137) 

1 

cm 

An  INPUT  parameter  used  to  calculate 

the  DY  array  if  the  axial  dimension 
of  the  cells  is  uniform. 

EAMMP 

(PH2) 

1 

ergs/g 

Specific  internal  energy  of  mass 

transported  across  the  right  edge 
of  cell. 

EAMPY 

(PH2) 

1 

ergs/g 

Specific  internal  energy  of  mass 
transported  across  top  of  cell. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

ECK 

■2(76} 

1 

-- 

Used  in  EDIT.  Relative  error  in 
energy  sum: 

■  ETHj/ETH; 

where  Ei.  is  total  energy  in  cell  K. 

If  lECKy  >  DMIN,  execution  is 
stopped.' 

EMIN 

“Z(8S} 

1 

ergs/g 

INPUT  parameter.  Minimum  specific 
internal  energy  of  an  ideal  gas  to 
be  used  in  the  pressure  iteration 
in  CDT.  Usually  EMIN  =10'. 

EMOB 

■ZCISO) 

1 

ergs 

Calculated  in  PH2 .  Printed  in  EDIT. 
Energy  transported  across  bottom  of 
mesh. 

EMOR 

=  2(135} 

1 

ergs 

Calculated  in  PH2.  Printed  in  EDIT. 
Energy  transported  across  right  side 
of  mesh. 

EMOT 

=2(146} 

1 

ergs 

Calculated  in  PH2.  Printed  in  EDIT. 
Energy  transported  across  top  of 
mesh. 

ENERGY 

B.C. 

1 

ergs/g 

Defined  in  CDT  as  specific  internal 
energy  of  a  material.  Used  in  EQST 
where  P  =  f (ENERGY, RHOW) . 

EOB 

=2(134} 

1 

ergs 

Calculated  in  PHI.  Printed  in  EDIT. 
Energy  change  due  to  work  done  at 
bottom  boundary  of  grid. 

EOR 

=2(132} 

1 

ergs 

Calculated  in  PHI,  Printed  in  EDIT. 
Energy  change  due  to  work  done  at 
right  boundary  of  grid. 

EOT 

-2(133} 

1 

ergs 

Calculated  in  PHI.  Printed  in  PHTT. 
Energy  change  due  to  work  done  at  top 
boundary  of  grid. 

ERDUMP 

B.C. 

1 

... 

Used  in  EDIT  and  ERROR.  Flags  EDIT 
to  do  only  a  tape  dump  on  an  error 
exit. 

ESA 

(EQST} 

30 

Defined  in  DATA  Statement.  Values 
of  "a"  in  equation  of  state  for  19 
materials.  [=(Y-l}when  using  perfect 
gas  equation  of  state.  ] 

ESALPH 

(EQST} 

30 

Defined  in  DATA  Statement.  Values  of 
"o"  in  equation  of  state  for  19 
materials. 

ESB 

(EQST} 

30 

Defined  in  DATA  Statement.  Values  of 
"b"  in  equation  of  state  for  19 
materials. 

ESBETA 

(EQST} 

30 

Defined  in  DATA  Statement.  Values  of 
"B"  in  equation  of  state  for  19  materials 

ESCAPA 

(EQST} 

30 

dynes/cm^ 

Defined  in  DATA  Statement.  Values  of 

A"  in  equation  of  state  for  19  materials 
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Variable  Number  of 

Name  Location  Words  Units 


Definition 


ESCAPfi 

CEQST] 

30 

2 

dynes/cm 

Defined  in  DATA  Statement.  Values  of 
"B"  in  equation, of  state  for  19  materials. 

ESES 

(EQST) 

30 

ergs/g 

Defined  in  DATA  Statement.  Values 
of  "E  "  in  equation  of  state  for  19 
materials. 

ESESP 

(EQST) 

30 

ergs/g 

Defined  in  DATA; Statement.  Values  of 
"E'"  in  equation  of  state  for  19  materials 

ESEZ 

(EQST) 

30 

ergs/g 

Defined  in  DATA  Statement.  Values 
of  "E  "  in  equation  of  state  for 

19  materials. 

ETH 

=Z(13) 

1 

ergs 

Theoretical  value  of  total  energy  in 
the  mesh.  Calculated  in  SETUP  ini¬ 
tially;  updated! in  PHI,  PH3 ,  PHZ  and 

REZONE .  i 

EVAPEN 

=Z(101) 

1 

ergs 

Calculated  in  PHZ  and  CDT.  Printed 
in  EDIT.  Sura  of  energy  lost  through 
"evaporation"  of  mass  left  in  cells 
due  to  round-off  error.  Adjusted  in 

CDT  when  "evaporating"  isolated  cells. 
Initialized  in  SETUP. 

EVAPM 

=Z(100) 

1 

g 

Calculated  in  PHZ  and  CDT.  Printed 
in  EDIT.  Sum  of  mass  lost  through 
"evaporation".  See  EVAPEN. 

EVAPMU 

=Z(102) 

1 

cm-g/sec 

Calculated  in  PHZ  and  CDT.  Printed 
in  EDIT.  Sum  of  radial  momenta  lost 
through  "evaporation'.'.  See  EVAPEN. 

EVAPMV 

=Z(1Q3) 

1 

cm-g/sec 

Calculated  in  PHZ  and  CDT.  Sum  of 
axial  momenta  lost  through  "evapora¬ 
tion".  See  EVAPEN. 

EZPH2 

=Z(104) 

1 

ergs/g 

Sum  of  specific  internal  energy  fluxes 
that  are  less  than  UMIN*  and  are  set 
to  zero  in  PHZ.  Printed  in  EDIT. 

FINAL  - 

=2(113) 

1 

INPUT  parameter.!  Maximum  value  of 
stability  fraction  (STAB) .If 

FINAL  •  0;  the  stability  fraction 
will  be  constant.  Used  in  CDT. 

FLEFT 

(PH2) 

100 

cro-g/sec 

Radial  momentum  of  mass  transported 
across  left  side  of  cell.  Equivalenced 
to  UL  array.  (See  Appendix  B) 

FRACRT 

B.C. 

1000 

cm 

The  fractional  area  of  the  right  cell 
boundary  associated  with  a  given 
material  in  a  mixed  cell.  Used  to 
compute  mass  flux  of  that  material. 

FRACTP  . 

B.C. 

1000 

cm 

The  fractional  area  of  the  top  cell_ 
boundary  associated  with  a  given 
material  in  a  mixed  cell.  Used  to 
compute  mass  flux  of  that  material. 

GAMC 

(PH2) 

100 

g 

Mass  transported  across  left  side  of 
cell.  Equivalenced  to  PL  array. 

(See  Appendix  B) 
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Variable  Ntunber  of 

Nane  Location  Words  Units 


Definition 


GAMMA 

■Z(62) 

1 

INPUT  parameter,  GAMMA  ■  ESA+1.  Used 
In  EQST  for  an  Ideal  gas.  User  should 
set  GAMMA  -  0.  when  problem  does  not 
Involve  an  ideal  gas. 

I 

•Z(88) 

1 

-- 

Used  in  most  subroutines  as  index  in 
radial  direction. 

ICSTOP 

■Z(7) 

1 

INPUT  parameter.  Used  in  EDIT.  Execu¬ 
tion  stops  on  ICSTOP  cycle  when  stop¬ 
ping  on  a  specified  cycle  rather  than 
time. 

GLUED 

»Z(79) 

1 

INPUT  parameter.  GLUED  =  1.  Causes 
all  mixed  cells  to  be  "glued"  Csee 
discussion  of  GLUE  subroutine  in  Sec¬ 
tion  7.1.10)  after  the  PHI  work  terms 
have  been  used  to  update  U,  V  and  E. 
This  option  is  rarely  used  but  can  be 
■  helpful  in  stabilizing  the  interface 
between  a  hot  gas  expanding  into  a  cold 
dense  material. 

lEXTX 

-Z(123) 

1 

Used  in  REZONE.  When  lEXTX  ■=  1.,  grid 
is  rezoned  in  x-direction.  Enables 
user  to  rezone  grid  in  x-direction 
only. 

IMAX 

'Z(33) 

1 

•* 

INPUT  parameter.  Number  of  columns 
in  mesh.  IMAX  must  be  an  even  number 
if  grid  is  to  be  rezoned  in  the  x- 
direction- 

IMAXA 

-Z(34) 

1 

— 

IMAX  +  1, 

INTER 

"Z(87) 

1 

INPUT  parameter.  On  print  cycles  if 
INTER  !<  0,  EDIT  will  print  after  CDT, 
PHI,  and  PH3,  If  INTER  =99,  in 
addition  to  extra  EDIT  prints,  energy 
total  and  elastic  and  plastic  work 
are  printed  in  PH3.  If  INTER  =  7, 
energy  totals  are  printed  in  PH2  in 
addition  to  the  extra  EDIT  prints. 

IPCYCL 

-Z(49) 

1 

INPUT  parameter.  Used  in  EDIT.  The 
number  of  cycles  between  EDIT  prints 
when  editing  on  cycles  rather  than 
time . 

IPR 

-zcis) 

1 

Maximum  number  of  iterations  to  be 
performed  by  CDT  to  achieve  pressure 
equilibration  between  materials  in 
mixed  cells. 

IVARDX 

=  Z(83) 

1 

-- 

When  IVARDX  =  1,  radial  dimension 
of  ctnis  is  variable. 

IVARDY 

-Z(54) 

1 

-- 

When  IVARDY  =  1,  the  axial  dimen¬ 
sion  of  cells  is  variable. 
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Variable 
.  Name 

Location 

Number  of 
Words 

Units 

Definition 

11 

-2(47) 

1 

-- 

INPUT 

limit 

parameter.  11  is  used  to 
calculation  in  the  radial 

direction  to  an  "active  mesh."  Be¬ 
yond  II  nothing  is  yet  distrubed  from 
the  initial  conditions.  II  is 
specified  initially  as  (2  +  the 
column-number  of  the  last  column 
in  which'  there!  is  a  non-zero  velo¬ 
city  or  internal  energy).  II  is 
increased  automatically  as  inactive 
cells  become  active.  However.  II  is 
never  larger  than  IMAX. 


12 

-Z(48) 

1 

•• 

INPUT  parameter.  Like  11  but  for 
axial  disturbance  limit. 

13 

B.C. 

1 

Used  in  EDIT  as  a  flag  for  "short" 
or  "long"  prints.  13  =  number  of 
columns  of  grid  whose  U,  V,  P,  AIX, 
etc.  are  listed. 

J 

-2(89) 

1 

-- 

Used  as  row-index  in  most  subroutines 

JCENTR 

=2(147) 

1 

Used  in  GLUE  and  PH2  to  maintain 
symmetry  of  calculation  involving  a 
sphere.  JCENTR  =  J-line  coincident 
with  sphere  center. 

JEXTY 

-2(124) 

1 

-- 

A  flag  used  in  REZONE.  When  JEXTY  = 
grid  is  rezoned  in  y-direction. 

JMAX 

=  2(3S) 

1 

-- 

INPUT  parameter.  Number  of  rows  in 
the  grid.  JMAX  must  be  an  even 
number  if  the  grid  is  to  be  rezoned. 

JMAXA 

to 

II 

1 

— 

Calculated  in  SETUP,  (JMAX  •►1). 

K 

0 

o 

1 

.  .. 

Used  as  cell-index  in  all  subroutines 
K  -  (J-1)»IMAX  +  I  +  1. 

KA 

B.C. 

1 

-- 

Used  as  cell-index  for  cell  above. 

KB 

B.C. 

1 

-- 

Used  as  cell-index  for  cell  below. 

KMAX 

-2(37) 

1 

-- 

Calculated  in  SETUP  (=  IMAX*JMAX  +  1) 
K-index  of  last  cell  in  grid. 

KMAXA 

-2(38) 

1 

-- 

Calculated  in  SETUP  (KMAX  +  1) . 

KSPACE 

(EDIT) 

1 

-- 

Used  for  spacing  printed  output. 

KUNITR 

-2(14) 

1 

•• 

Unit  number  of  tape  read  by  INPUT. 
Must  be  7  when  reading  a  "CLAM" 
tape  (i.e.,  when  PK(3)  =  -3.). 

KUNITW 

-2(17) 

1 

Unit  number  of  tape  written  on  by 
SETUP  and  EDIT.!  Set~RWm  i<  7 
when  reading  a  "CLAM"  tape  to  pre¬ 
vent  SETUP  from,  writing  over  it. 

M 

■Z(91) 

1 

-- 

Usually  used  as  an  index  for  cell 

K  in  the  mixed  cell  arrays.  M  - 
MFLAG(K)-100. 

118 


3SR-350 


Variable  Number  of 

Name  Location  Words  Units  Definition 


MA 

B.C. 

1 

Usually  value  of  MFLAG  for  cell  above 
cell  K. 

MAPS 

-Z(42) 

1 

INPUT  parameter  that  determines  the 
printing  of  symbolic  maps  on  EDIT 
cycles.  If  MAPS  =0,  maps  are  not 
printed;  if  MAPS  =  1,  they  are 
printed.. 

MAT 

B,C. 

30 

Material  code  numbers  for  the  packages 
MAT(l)  =  3  indicates  the  material  in 
package  one  is  iron.  See  list  in 
Section  2.3  of  this  report  or  com¬ 
ments  in  listing  of  EQST  subroutine. 

MAXX 

=Z(120) 

1 

— 

Used  in  ADDTCR.  Right  column  and  top 
row,  respectively,  of  region  in  which 

MAXY 

•  -Z(122) 

1 

■  ** 

tracer  particles  are  to  be  added 
(see  NADD,  MINX,  MINY] . 

MBS 

»Z(116)  . 

1 

Maximum  number  of  coordinates  used  to 
describe  a  package  which  extends 
beyond  the  original  grid.  (Maximum 
value  in  MPAC  array] . 

MBBB 

-zcns) 

1 

•• 

Number  of  material  packages  set  up 
with  finite  boundaries  extending  be¬ 
yond  the  limits  of  the  original  grid. 

MFK 

B.C. 

1 

” 

Usually  =  MFLAG (K]. 

MFLAG 

B.C. 

2500 

Flag  for  each  cell  that  indicates 
whether  it  is  pure  or  mixed.  If 
MFLAG(K]  <  100  cell  is  pure.  If 
MFLAG(K)  >  100  cell  is  mixed  and 
flag  gives  storage  location  in  mixed 

cell  arrays  for  cell  K  (see  SIE). 

MINX 

•Z(119) 

1 

” 

[Used  in  ADDTCR.  Left  column  and 

1  bottom  row  resnectively  of  region 

MINY  . 

-Z(121) 

1 

in  which  tracer  particles  will  be 
[added  (see  NADD,  MAXX,  MAXY]. 

MO 

B.C. 

1 

Defined  in  INFACE.  MFLAG (K]  of 
cell  K  before  it  becomes  mixed. 

Used  in  NEWMIX  to  define  mixed  cell 
variables  for  cell  K. 

MPAC 

B.C. 

3 

MPAC(N]  is  the  number  of  coordinates 
used  to  describe  the  boundaries  of 
material  package  N  which  extends 
beyond  the  boundaries  of  the  original 
grid. 

MPACK 

B.C. 

3 

•  • 

MPACK(M)  is  the  package  number  of  the 
Mth  package  that  extends  beyond  the 
boundaries  of  the  original  grid.- 

MR 

B.C. 

1 

— 

Value  of  MFLAG  for  cell  on  right. 

N 

B.C. 

1 

-- 

In  EQST,  N  is  material  code  number 
transferred  from  CDT. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

N6 

»Z(56) 

1 

INPUT  parameter.  Used  in  EQST. 
Negative  pressures  are  allowed  in 
cells  above  J  =  N6.  N6  is  redefined 
in  RfZONE.  Setting  N6  »  o  allows 
negative  pressures  everywhere.  On 
the  other  hand  to  make  sure  that 
negative  pressures  are  always  set 
to  zero  give  N6  a  very  large  value 
(many  times  as  large  as  JMAX] . 

NIO 

-Z(60) 

1 

-- 

Used  in  CDT  to  identify  I- index  of 
cell  which  controls  the  time  step. 

Nil 

-Z(61) 

1 

-- 

Used  in  CDT  to  identify  J-index  of 
cell  which  controls  the  time  step. 

NADD 

=Z(118) 

1 

-- 

INPUT  parameter.  ■  When  NADD  >  0  on 

the  first  cycle  of  a  run,  ADDTCR 
is  called  and  adds  at  most  NADD 
material  tracer  particles  between 
existing  tracers  in  a  given  region. 
It  will  not  put  more  than  NTRACR 
tracers  in  a  cell  for  each  material 
(see  also  MINX,  MINY,  MAXX,  MAXY) . 


NC 

=  Z(30) 

1 

Cycle  number.  Set  initially  to  -1 
in  input.  Incremented  thereafter 
in  CDT. 

NDUMP7 

-Z(6] 

1 

INPUT  parameter.  ■  Used  in  EDIT  to 
control  frequency  of  tape  dumps. 

E.g.,  tape  dump  will  occur  on  every 
EDIT  print  when  NDUMP7  =  1,  or  on 
every  fifth  EDIT  print  when  NDUMP7 
=  5. 

NECYCL 

=  Z(77) 

1 

-• 

Defined  and  printed  in  EDIT.  The 
cycle  during  which  the  largest 
relative  error  in  the  energy  sum 
was  computed. 

NERR 

B.C. 

1 

Used  in  ERROR  as  exit  flag.  Prevents 
ERROR  from  being  called  more  than 
once  during  a  single  run. 

NFRELP 

■2(5) 

1 

INPUT  parameter.  Used  in  EDIT  to 
control  frequency  of  "long"  prints. 

A  "long"  print  will  occur  with  every 
EDIT  if  NFRELP  =  1;  with  every  fifth 
EDIT  if  NFRELP  =  5  (see  13). 

NK 

B.C. 

1 

Tells  which  statement  number  of  a 
subroutine  is  near  the  call  to 

ERROR.  i 

NMAT 

-Z(68) 

1 

INPUT  parameter.  'Equals  number  of 
material  packages'  in  the  problem. 
[Note;  each  material  package  can 
be  made  up  of  several  disconnected 
subpackages.  Also  package  N  and 

N-*-!  can  contain  the  same  material 

but  Have  different  initial  conditions 
(P.V.E)].  I 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

NMXCLS 

-ZC73) 

1 

•  • 

INPUT  parameter.  Equals  the  maxi¬ 
mum  number  of  mixed  cells  that 

SETUP  or  NEWMIX  can  generate.  This 
number  should  coincide  with  the 
dimension. of  variables  in  the 

MXCELL  common  block. 

NODUMP 

■Z(96) 

1 

INPUT  parameter.  Used  in  EDIT. 

When  NODUMP  =  1 ,  no  tape  dumps  are 
made  except  by  SETUP  on  cycle  0. 
Allows  user  to  restart  a  problem 
without  writing  on  the  restart  tape. 

NPRINT 

B.C. 

.  1 

NPRINT  "  1  during  the  cycle  on  which 
EDIT  prints  and  checks  the  energy 
error. 

NR 

B.C. 

1 

“  “ 

Identifies  which  subroutine  calls 
ERROR.  Used  in  ERROR  for  printing 
error  message. 

NRC 

•Z(3Z] 

1 

-- 

Used  in  PHI  and  PH2  in  advancing 
active  grid. 

NREZ 

-Z(20) 

1 

-- 

Defined  in  SETUP.  Equals  maximum 
number  of  rezones  allowed. 

NTCC 

-Z(81] 

1 

•• 

INPUT  parameter.  When  NTCC  =  1 
passive  cell  centered  tracers 
are  used. 

NTRACR 

-Z(72) 

1 

INPUT  parameter.  Used  to  determine 
density  of  material  tracer  particles, 
If  NTRACR  »  2,  SETUP  will  place  two 
tracers  per  cell  edge  around  a 
material  package?. 

NS IDES 

-Z(67) 

1  . 

INPUT  parameter.  Equals  number  of 
line  segments  —  other  than  grid 
lines  —  that  bound  the  void.  If 
there  is  no  void  in  the  problem, 
NSIDES  =  0.  Used  in  SETUP  to  read 
in  cards  for  defining  the  tracer 
particles  that  circumscribe  the 
void  package. 

NTPMX 

-ZC78) 

1  ■ 

INPUT  parameter.  Equals  the  maxi¬ 
mum  number  of  tracers  that  SETUP 
or  ADDTCR  can  generate  to  circum¬ 
scribe  a  single  material  package. 

The  value  of  this  number  should 
coincide  with  the  dimension  of  the 

TX  and  TY  arrays . 

NUMREZ 

-Z(12) 

1 

INPUT  parameter.  Initially  equals 
number  of  rezones  allowed  in  one  . 
run.  Diminished  by  one  after  each 
rezone.  Used  in  EDIT  and  REZONE. 

NUMSCA 

-Z(43) 

1 

-- 

INPUT  parameter.  Number  of  times 
the  print  interval  is  to  be  re¬ 
scaled.  Used  in  EDIT.  See  PRDELT 

for  further  details. 
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Variable 

Name 

NUMSP 

NVOID 

NZ 

P 

PACX 

PACY 

PL 

PLW 

PMIN 

PRCNT 

PRDELT 


Number  of' 

Location  Words  Units 


Definition 


=2(4) 


B.C. 


-ZC19) 


B.C. 


B.C. 


B.C. 


B.C. 


MXCELL 


-ZC86) 


=ZC16) 


■2(45) 


1 


1 


1 


2500  dynes/cm^ 

30  cell 

30  cell 

2 

200  dynes/cm 

3  ergs 

1  dynes/cm 


1 


1  sec 


Used  in  EDIT  to  count  the  number  of 
prints  (short  or  long)  since  the 
last  tape  dump. 

Defined  in  INPUT:  =  NMAT  +1, 
number  of  material  packages  ♦  1. 

The  number  of  the  void  package  if 
there  is  one. 

Defined  and  used  in  EDIT  for  1-D 
problems.  NZ  =  4**NRZ.  After  re¬ 
zoning  the  grid  NZ  is  used  to  scale 
the  values  printed  by  EDIT  for  the 
total  mass,  energy  and  momentum. 

Cell-pressure.  ,  (IMAX  by  JMAX  array). 
Calculated  in  EQST.  Used  by  PHI. 
[Note:  the  P-storage  space  is  used 

tor  other  information  in  INFACE,  PH3 
and  PH2.] 

PACXCN,M)  is  the  x-coordinate  of  the 
Mth  tracer  of  the  Nth  package  which 
has  boundaries  extending  beyond  the 
original  grid. 

PACY(N,M)  is  the  y-coordinate  of  the 
Mth  tracer  of  the  Nth  package  which 
has  boundaries  extending  beyond  the 
original  grid.  ' 

Used  in  PHI  for,  saving  pressures  on 
left  side  of  cell.  (Storage  used  for 
other  purposes  in  CDT  and  PH2.) 

An  array  that  stores  the  total  work 
done  by  plastic;  stresses  in  each 
material  package.  Calculated  in  PH3; 
printed  in  EDIT. 

Used  as  a  pressure  cut-off.  Cal¬ 
culated  and  printed  in  CDT:  PMIN  = 
min  (PMIN,  C'  p,  UMIN) .  For  cycle  0, 
PMIN  defined'by  an  input  card. 

Convergence  requirement  for  equiliV 
brating  pressures  in  a  mixed  cell. 

/•’i  • 

If  { - — —  W  PRCNT  for  all  materials 

(i)  in  cell  K,  materials  in  cell  are 
considered  to  be  pressure  equilibrated 
and  P(K)  -  7. 

INPUT  parameter.  Gives  the  initial 
time  interval  between  EDIT  prints. 
There  are  five  parameters  which  con¬ 
trol  printing  frequency:  PRDELT, 
IPCYCL,  PRLIM,  PRFACT,  and  NUMSCA. 

If  the  user  is  printing  on  time 
(PRDELT  )<  0  and  IPCYCL  =  0)  ,  DT  will 
be  adjusted  so  that  a  print  will 
occur  exactly  every  PRDELT  seconds. 

If  the  user  is  printing  on  cycles 
(PRDELT  =  0,  IPCYCL  )<  0)  ,  a  print 
will  occur  every  IPCYCL  cycles. 


122 


3SR-350 


Variable  Number  of 

Name  Location  Nords  Units 


PRDELT 

(continued] 


PRESUR 

B.C.  . 

1 

dynes/cm 

PRFACT 

-Z(46] 

1 

-- 

PRLIM 

-Z(44] 

1- 

— 

PROB 

=  Z(1] 

1 

— 

PROP 

B.C, 

50 

— 

PRR 

(PHI) 

1 

dynes/cm' 

PRTIME 

-Z(131] 

1 

sec 

RADIUS 

-2(148) 

1 

cm 

RATIO 

(CDT] 

1 

sec 

RC 

(PHI) 

1 

cm 

RELERR 

(EDIT) 

1 

Definition 


PRLIM,  PRFACT  and  NUMSCA  are  used  to 
increase  the  print  interval.  PRLIM 
is  the  time  (or  cycle]  at  which  PRDBLT 
(or  IPCYCL]  and  PRLIM  are  multiplied 
by  PRFACT.  The  new  value  of  PRLIM 
establishes  the  next  tine  (or  cycle) 
when  the  print  interval  will  be  re¬ 
scaled.  This  process  continues  at  most 
NUMSCA  times . 

EXAMPLE:  You  wish  to  print  every 
1  10'  sec  until  you_yeach  1  x  10' 

sec,  thfn  every  1  x  lo"  sec  ugtil 
1  *  10"  sec  and  every  1  x  lo'  sec 
thereafter: 

PRDELT  -  1.  X  10'; 

PRLIM  =  1.  *  10' 

PRFACT  »  10. 

NUMSCA  ■=  2. 

Defined  in  EQST:  pressure  =  f(p,E]. 
Used  in  CDT  to  define  P(K)  in  the  case 
of  pure  cells,  and  in  the  case  of  mixed 
cells  to  define  the  pressure  of  a 
material . 

INPUT  parameter.  Used  in  EDIT  as  a 
factor  for  rescaling  PRDELT,  IPCYCL, 
and  PRLIM  when  PRLIM- time  or  cycle  is 
reached  (see  PRDELT].  Should  be  >  1. 

INPUT  parameter:  time  or  cycle  at 
which  to  rescale  PRDELT  (or  IPCYCL) 
and  PRLIM  by  PRFACT  (see  PRDELT). 

INPUT  parameter.  .  Identifying  problem 
number.  Should  be  between  0  and  100. 

Calculated  and  used  in  EDIT.  For  1-E 
problems  the  totals  for  energy,  mass 
momentum  per  unit  area  are  computed 
by  dividing  by  4**  (number  of  rezones) . 
These  totals  are  stored  in  PROP  for 
printing.  Equivalenced  to  UL  array. 

Pressure  at  right  cell  boundary. 

Initially  set  to  PRDELT  in  INPUT. 
Thereafter  calculated  in  EDIT.  When 
T  •  PRTIME,  it  is  time  to  print. 

INPUT  parameter.  Radius  of  sphere 
(in  cms) . 

Used  in  calculation  of  DT:  ratio  of 
speedljji^  for  a  given  cell. 

Distance  from  axis  to  center  of  a  cell. 

Used  for  storing  and  printing  maximum 
relative  error  in  the  energy  sum. 
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Variable  Number  of 

Name  Location  Words 


Units 


Definition 


REZ  =Z(95)  1 


REZFCT 

=  Z(71) 

1 

•  - 

RHOIN 

B.C. 

3 

g/cm^ 

RllOW 

B.C. 

1 

g/cm^ 

RHOZ 

MXCELL 

30 

g/cm^ 

RMU 

B.C. 

3 

dynes/cm' 

ROEPS 

=Z(110) 

1 

-- 

RR 

(PHI) 

1 

cm 

RTM 

=Z(57) 

1 

g 

RTMU 

=2(10) 

1 

cm-g/sec 

RTMV 

=2(58) 

1 

cm-g/sec 

SAMMP 

MXCELL 

750 

g 

SAMPY 

MXCELL 

750 

g 

Flag  defined  in  PH2  and  used  in  EDIT. 
Signals  when  the  REZONE  subroutine 
should  he  called.  'The  rezone  flag  is 
turned  on  when  material  begins  to  flow 
out  through  transmittive  boundaries 
(see  VT) .  If  the  grid  is  rezoned  in 
both  directions,  each  set  of  four  cells 
in  the  mesh  is  made  into  one  cell.  The 
new  mass  is  the  sum  of  masses  in  the  four 
original  cells.  Momentum  and  total 
energy  are  conserved  but  in  so  doing 
some  kinetic  energy  is  changed  to  in¬ 
ternal.  (The  result  is  that  rezoning 
has  a  stabilizing  influence.)  When 
all  rezones  have  been  performed,  or 
if  no  rezones  are  permitted,  material 
is  allowed  to  flow  out  through  trans¬ 
mittive  boundaries  and  the  mass  and 
energy  are  lost  from  the  system. 

INPUT  flag  for  rezoning.  If  =  1.,  the 
grid  is  rezoned  (NUMREZ)  times.  If  =  0., 
no  rezoning  is  done.  Tested  in  PH2  and 
EDIT. 

The  initial  densities  of  material  in 
the  packages.  Defined  by  input  cards 
read  by  SETUP. 

Density  of  material.  Defined  in  CDT 
and  used  in  EQST  to  define  pressure: 

P  =  f (ENERGY, RMOW) . 

Defined  in  DATA  statement  in  BLOCK. 

Normal  density  for  'IP  materials. 

Rigidity  modulus  values  for  material 
packages.  Defined  by  input  cards  read 
by  SETUP.  Used  in  PII3. 

INPUT  parameter.  A  round-off  epsilon 
used  in  calculating  cutoffs  of  energy, 
velocity,  pressure  and  mass  flux. 

Distance  from  axis  to  center  of  cell 
on  the  right. 

Calculated  in  PH2.  Printed  in  EDIT. 
Total  mass  lost  out  right  side  of  grid. 

Calculated  in  P1I2.  Printed  in  EDIT. 
Total  radial -momentum  lost  out  right 
side  of  grid. 

Calculated  in  PU2.  '  Printed  in  EDIT. 
Total  axial-momentum  lost  out  right 
side  of  grid. 

The  flux  of  materials  across  right 
boundary  of  mixed  cells.  Calculated 
in  FLUX;  used  in  Pli2  and  INFACE. 

The  flux  of  materials  across  the  top 
boundary  of  mixed  cells.  Calculated  in 
FLUX;  used  in  P112  and  INFACE. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

SDT 

B.C. 

1 

sec 

Defined  in  INFACE  »  DT/CYCMX.  Time  step 
for  each  subcycle  of  INFACE. 

SIE 

MXCELL 

7S0 

ergs/g 

The  specific  internal  energy  of  materials 
in  mixed  cells.  If  cell  K  is  mixed  and 

M  •  MFLAG{K)-100,  then  SIE(1,M)  is  the 
specific  internal  energy  of  material  one 
in  cell  K  (see  XMASS,  RHO) . 

SIENEW 

(PH2) 

■  1 

ergs/g 

New  value  of  specific  internal  energy. 

SIGC 

B.C. 

100 

ergs 

Used  in  PH2  for  internal  energy  of  mass 
transported  across  left  side  of  cell. 
Equivalenced  to  PL (103).  (See  Appendix  B) 

SIGMU 

(PH23 

1 

cm-g/ sec 

Total  change  in  radial  momentum  of 
a  cell. 

SIGMV 

(PH2) 

1 

cra-g/sec 

Total  change  in  axial  momentum  of  a 
cell. 

SNB 

B.C. 

so 

dynes/cm^ 

Normal  deviator  stress  at  bottom  of  a 
cell.  Equivalenced  to  P  array. 

SNL 

(PH3) 

1 

dynes /cm ^ 

Normal  deviator  stress  at  left  of  a  cell. 

SNLX 

CPH3) 

1 

dynes /cm  . 

[-  SNL»X(I-1)]. 

SNR 

(PH3) 

1 

dynes/cm^ 

Normal  deviator  stress  at  right  of  a 
cell.  (“Sjj  at  right). 

SNT 

(PH3} 

1 

dynes/cm^ 

Normal  deviator  stress  at  top  of  a 
cell  ('822  f°P) • 

SOLID 

(STRNG) 

1 

g/cm^ 

=  RHOZ^  »  AMDMi.  In  PH3,  if  PUOi  <  SOLID 
for  one  material,  i,  in  cell  K,  stresses 
of  cell  K  are  set  to  zero. 

SRATIO 

CCDT) 

1 

sec 

Used  to  calculate  DT.  The  smallest  ratio 
in  grid  of  a  cell's  minimum  dimension  to 
the  sum  of  its  maximum  velocity  and  sound 
speed. 

SSI 

-  2(127) 

1 

g/ergs 

Used  in  EQST:  in  the  expanded  equation 
of  state:  =  1/(EJ  -  Eg) . 

SS2 

=  2(128) 

1 

Constant  controlling  reflective  (and 
axis)  boundary  treatment.  Defined  in  . 

INPUT;  SS2  =  1.  Used  in  PH2. 

SS4 

-  2(130) 

1 

INPUT  parameter.  If  SS4  f  0.,  REZONE  is 
called  on  first  cycle  of  a  run.  Used  in 
EDIT. 

SSIEN 

B.C. 

3 

ergs/g 

Initial  specific  internal  energy  of  each 
package.  Defined  by  input  cards  read  in 
SETUP. 

STAB 

-  2(139) 

1 

INPUT  parameter.  Used  in  CDT.  Initial 
value  of  "stability  fraction"  for  the 
calculation  of  DT.  If  FINAL  »  0.,  STAB 
is  constant.  Otherwise  its  value  progresses 
from  STAB  to  FINAL  in  a  geometric  progfessioi 
[Note:  DT  •  STAB  •  SRATIO] 
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Variable  Number  of 


Name 

Location 

Words 

Units 

Definition 

STB 

(PH3) 

1 

2 

dynes/cn 

Shear  stress  at  bottom  of  cell. 

STEZ 

B.C. 

3 

ergs/g 

Eq  for  each  material  package.  Used 
in  yield-strength  calculation  in  STRNG. 
Defined  by  input  cards  read  by  SETUP. 
(See  STRENG) 

1 

STKl 

B.C. 

3 

dynes/ cm^ 

Y,  for  each  material  package.  Used 
in  yield-strength  calculation  in  STRNG. 
Defined  by  input  cards  read  by  SETUP. 
(See  STRENG)  1 

STK2 

B.C. 

3 

2 

dynes/cm 

Y,  for  each  material  package.  Used  in 
yield-strength  calculation  in  STRNG. 
Defined  by  input  cards  read  by  SETUP. 
(See  STRENG) 

STL 

(PH3} 

1 

dynes/cm^ 

Shear  stress  at  left  of  cell. 

STLX 

(PH3) 

1 

dynes/cm 

-  STL  *  X(I-l) . 

STR 

CPH3) 

1 

dynes/cm^ 

Shear  stress  at  right  of  cell. 

STRENG 

(STRNG) 

1 

2 

dynes/cm 

Yield  strength  of  material: 

STRENG  =CY,  +  Y,ii  +:YjU*).  (1  -  E/E,). 

U  *  p/p,  -  1.  If  STRENG  £  0.,  stresses 
are  set  to  0.  If  E  >  E,,  STRENG  =  0  and 
stresses  are  set  to;  0. 

Y  is  CZERO, 

0 

Y  is  STKl, 

1 

is  STK2, 

E  is  specific  internal  energy  of  material, 
E  is  STEZ, 

0  i 

p  is  density  of  material, 


is  normal  density  of  material. 


STRSRR 

ELPL 

2500 

dynes/cm^ 

The  cell  centered  deviatoric  normal  stress 
in  the  radial  direction  •  (IMAX  by  JMAX  array) 

STRSRZ 

ELPL 

2500 

2 

dynes /cm 

The  cell  centered  deviatoric  shear  stress 
in  the  radial  direction  (IMAX  by  JMAX  array) 

STRSZZ 

ELPL 

2500 

2 

dynes /cm 

The  cell  centered  deviatoric  normal  stress 
in  the  axial  direction  (IMAX  by  JMAX  array). 

STT 

(PH3) 

1 

2 

dynes/cn 

Shear  stress  at  top  of  cell. 

SUM 

B.C. 

1 

-- 

Used  in  most  routines  as  working  storage 
when  summing  quantities. 

SUME 

(PH2) 

1 

ergs 

Used  in  PH2.  Sums  energy  fluxes  ignored 
during  one  cycle.  Used  to  adjust  ETH, 
the  theoretical  energy  total. 

T 

•  Z(84) 

1 

sec. 

Time.  Initially  defined  in  INPUT.  Incre¬ 
mented  in  CDT.  Adjusted  in  EDIT  for 
printing  and  stopping  at  specified  time 
values. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

TAU 

B.C. 

SO 

cm^ 

Initially  defined  in  SETUP  and  redefined 
in  RE20NE.  Area  of  cell  face: 

»  jr{X(I)*  -  X(I-l)®].  Used  inmost 
subroutines. 

TAODTS 

(PHI) 

1 

ca?/sec 

-  TAU  *  DT. 

TH03 

CPH3) 

1 

.  1/sec 

-  1/3  (u^  *  V,  *  i). 

TNOW 

(EDIT) 

1 

sec 

Used  in  EDIT:  =  time  now;  saved  when 

EDIT  calls  CDT  after  calling  RE20NE. 

TOPM 

.  Z(63) 

i 

g 

Calculated  in  PH2.  Printed  in  EDIT. 

Total  mass  transported  across  top  of 
grid. 

TOPMU 

-  Z(9) 

1 

cm-g/sec 

Calculated  in  PH2.  Printed  in  EDIT. 

Total  radial -momentum  transported  across 
top  of  grid. 

TOPMV 

-  Z(66) 

1 

cm-g/sec 

Calculated  in  PH2.  Printed  in  EDIT. 

Total  axial -momentum  transported  across 
top  of  grid. 

TRIAL,  , 

Local 

1 

cm/sec 

Used  in  CDT.  Maximum  sum  in  the  entire 
grid  of  a  cell's  sound-speed  and  largest 
velocity  component.  Printed  in  CDT  as 
MAXCUV. 

TSTOP 

•  2(50) 

1 

sec 

INPUT  parameter.  Value  of  T  at  which 
execution  stops  when  stopping  on  time 
rather  than  cycles. 

TWOPI 

B.C. 

1  • 

- -  .. 

Calculated  in  INPUT.  Used  in  PH3: 

“  2ir. 

TX 

B.C. 

1600 

cell 

X-coordinates  of  tracer  particles  cir¬ 
cumscribing  the  material  packages. 

TY 

B.C. 

1600 

cell 

Y-coordinates  of  tracer  particles  cir¬ 
cumscribing  the  material  packages. 

U 

B.C. 

2500 

cm/sec 

Radial  velocity  of  cell  (IMAX  by  JMAX 
array) . 

uamMp 

(PH2) 

1 

cm/sec 

Radial  velocity  of  mass  transported  across 
right  cell  edge. 

UAMPY 

(PH2) 

1 

cm/sec 

Radial  velocity  of  mass  transported  across 
top  cell  edge. 

UEFF 

(INFACE) 

1 

cm/sec 

Radial  velocity  of  a  tracer  particle, 
computed  from  an  area-weighted  averaging 
scheme. 

UK 

B.C. 

150 

cm/ sec 

Temporary  storage  for  three  rows  of 
radial  velocities. 

UKT 

(PH3) 

1 

cm/sec 

Temporary  storage  for  U(K)  when  computing 
elastic-plastic  work  done  on  cell  'K. 

VL 

(PHI) 

200 

cm/sec 

Radial  velocity  at  left  boundary  of  cell 

K.  UL  array  equivalenced  to  the  FLEFT 
and  YAMC  arrays  used  in  PH2. 

UMIN 

-  2(129) 

1 

cm/sec 

Calculated  in  CDT.  Used  as  velocity  cutoff 
in  PH2  •  (ROEPS)  »  (sound-speed  + 
velocity)„„. 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Description 

UNxxx 

=  Z(xxx) 

1 

Unused  Z-storage. 

UNEW 

(PH2) 

1 

cm/sec 

New  radial  velocity  of  cell  K. 

uox 

(PH3]  . 

1 

1/sec 

”2U„/(Xi  ♦  Xj.i).  Used  to  define 
strain  rates  of  cell  K. 

URR 

B.C. 

1 

cm/sec 

Used  in  PHI,  FLUX  and  PH2.  Temporary 
storage  for  radial  transport  velocity 
at  right  boundary  of  cell  K. 

UUR 

B.C. 

3 

cm/sec 

Initial  radial  velocities  of  packages. 
Defined  by  input  cards  read  in  SETUP. 

UVMAX 

-  Z(80) 

1 

cm/sec 

Calculated  in  COT,  Maximum  velocity 
in  grid.  Redefined  in  PHI: 

/  \ 

max(uVMAX..^/E— j. 

Used  to  flag  cells  that  should  be  glued 
after  PHI. 

V 

B.C. 

2500 

cm/sec 

Axial  velocity  of  cell  (IMAX  by  JMAX 
array) . 

VABOVE 

B.C. 

1 

cm/sec 

Used  in  PHI,  FLUX,  and  PH2 ,  Temporary 
storage  for  axial  transport  velocity  at 
top  boundary  of  cell. 

VAMMP 

(PH2) 

1 

cm/sec 

Axial  velocity  of  mass  transported 
across  right. cell  edge. 

VAMPY 

(PH2) 

1 

cm/sec 

Axial  velocity  of  mass  moving  across 
top  cell  edge. 

VBLO 

(PHI) 

1 

cm/sec 

Axial  velocity  at  bottom  boundary  of 
cell. 

VEFF 

(INFACE) 

1 

cm/ sec 

Axial  velocity  of  a  tracer  particle, 
computed  from  an  a,rea-weighted  averaging 
scheme. 

VEL 

(PHI) 

1 

-- 

Subcycle  flag:  VEL  =  1.  on  first  pass 
through  PHI,  and  U,  V,  E  are  updated. 

On  second  pass  VEL  =  0.  and  only  E  is 
updated  using  new  U,V  from  first  pass. 

VK 

B.C. 

ISO 

cm/sec 

Temporary  storage  for  three  rows  of 
axial  velocities. 

VKT 

(PH3) 

1 

cm/ sec 

Temporary  storage  for  V(X)  when  computing 
elastic-plastic  work  done  on  cell  K. 

VNEW 

(PH2) 

1 

cra/sec 

New  axial  velocity  of  cell  K. 

VOLSPH 

(SETUP) 

1 

cm^ 

Volume  of  that  part  of  the  toroid  - 
generated  by  a  cell  which, is  inside  the 
sphere. 

VOW 

(EQST) 

1 

-- 
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Variable 

Name 

Location 

Number  of 
Words 

Units 

Definition 

VT 

-  Z(55) 

1 

g/cm^ 

Used  in  PH2  to  trigger  a  retone  of  the 
grid.  If  mass  is  transported  across 
grid  boundary  and  if  minimum  density  of 
mass  flux  (when  spread  over  a  cell) 

>  VT,  then  rezone  flag  (REZ)  is  set 

To  one.  Defined  in  INPUT:  VT  =  10'®. 

WA 

B.C. 

3; 

cm/sec 

Initial  axial  velocities  of  packages. 
Defined  by  input  cards  read  by  SETUP. 

WFLAGF 

■  ZCSl} 

1 

Used  in  INPUT  and  EDIT,  Set  =  1.  on 
first  cycle  of  run  in  INPUT.  Triggers 
an  EDIT  print  bn  first  cycle  of  every 
run.  Reset  to  0.  at  end  of  EDIT. 

WFLAGL 

=  Z(S2) 

1 

.  r  “ 

Used  in  MAIN  and  EDIT.  Flags  last  cycle. 
Set  '  1.  in  EDIT.  Signals  HELP  to  call 
exit. 

WS 

B.C. 

WSA  .  .  ■ 

B.C. 

WSB 

WSC 

B.C. 

B.C. 

1 

-- 

Used  in  most  subroutines  for  local 
working  storage. 

WSX 

B.C. 

WSY 

B.C. 

X 

B.C. 

SO 

cm 

Distance  from  axis  to  outside  of  cell. 
Equivalenced  to  XX  array  such  that  X(0) 

-  XX(1). 

XIENRG 

■ 

»-• 

o 

1 

ergs 

Total  internal  energy.  Calculated  in 

EDIT  and  used  for  printing  labels  on 
tracer  particle  plots. 

XKENRG 

=  Z(141) 

1 

ergs 

Total  kinetic  energy.  Calculated  in 

EDIT  and  used  for  orintina  labels  on 
tracer  particle  plots. 

XMASS 

MXCELL 

7S0 

g 

The  mass  of  materials  in  mixed  cells 
(see  SIE,  RHO) . 

XMAX 

■  Z(18) 

1 

cm 

-  X(IMAX). 

XP 

TRACRS 

1000 

cell 

X-coordinates  of  passive  cell-centered 
tracer  particles. 

XTENRG 

-  Z(142) 

1 

ergs 

Total  energy.  Calculated  in  EDIT  and 
used  for  printing  labels  on  tracer  particle 
plots . 

XUM 

(MAP) 

41 

•  • 

Used  in  MAP.  Array  has  negative  alpha¬ 
betic  characters  for  maps.  Defined  in 
a  DATA  statement. 

XX 

B.C. 

S2 

cm 

Equivalenced  to  X  array  so  as  to  make 

X(0)  available . 

Y 

B.C. 

100 

cm 

Distance  from  bottom  of  grid  to  top  of 
cell.  Equivalenced  to  YY  array  such 
that  Y(0)  -  YY(1). 
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Variable  Number  of 


Name 

Location 

Words 

Units 

Definition 

YAMC 

B.C. 

100 

cm-g/sec 

Calculated  and  used  in  PH2.  Axial 
momentum  of  mass  transported  across 
left  side  of  celli  Equivalenced  to  UL 
array.  (See  Appendix  B) 

YMAX 

(SETUP) 

1 

cm 

■=  Y(JMAX). 

YP 

TRACKS 

1000 

cell 

Y-coordinates  of, passive  cell-centered 
tracer  particles. 

YY 

B.C. 

102 

cm 

Equivalenced  to  Y  array  so  as  to  make 
Y(0)  available. 

Z 

B.C. 

1 

Storage  for  most, of  the  input  parameters 
follows  Z,  the  first  work  of  blank  commoi 
The  "Z-array"  (the  first  150  words  of 
blank  common)  is  written  on  tape  for 
restarting  problems  (see  Appendix  A). 
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APPENDIX  A 

Z-BLOCK  VARIABLES  LISTED  NUMERICALLY 


See  Dictionary  for  meaning  and  use. 


*1. 

PROB 

23. 

UN2  3 

*45. 

PRDELT 

2. 

CYCLE 

*24. 

DMIN 

*46. 

PRFACT 

3. 

DT 

25. 

UN25 

2*47. 

11 

4. 

NUMSP 

26. 

DTNA 

2*48. 

12 

2*5. 

NFRELP 

*27. 

CVIS 

2*49. 

IPCYCL 

2*6. 

NDUMP7 

28. 

UN28 

1*R50. 

TSTOP 

2*7. 

ICSTOP 

29. 

UN29 

51. 

WFLAGF 

8. 

PIDY 

30. 

NC 

52. 

WFLAGL 

9. 

TOPMU 

•  31. 

UN31 

53. 

UN5  3 

10. 

RTMU 

32. 

NRC 

2*54. 

IVARDY 

11. 

UNll 

2*33. 

IMAX 

55. 

VT 

2*12. 

NUMREZ 

34. 

IMAXA 

2*56. 

N6 

13. 

ETH 

2*35. 

UMAX 

57. 

RTM 

2*14. 

KUNITR 

36. 

JMAXA 

58. 

RTMV 

2*15. 

IPR 

37. 

KMAX 

59. 

UN59 

*16. 

PRCNT 

38. 

KMAXA 

60. 

NIO 

2*17. 

KUNITW 

39. 

BOTM 

61. 

Nil 

18. 

XMAX 

40. 

BOTMV 

*  62. 

GAMMA 

19. 

NZ 

41. 

NUMSPT 

63. 

■  TOPM 

20. 

NREZ 

2*42. 

MAPS 

64. 

BOTMU 

21. 

IGM 

2*43. 

NUMSCA 

65. 

UN65 

22. 

UN22 

*44. 

PRLIM 

66. 

TOPMV 
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2*67. 

NSIDES 

93. 

UN93 

2*119. 

MINX 

2*68. 

NMAT 

94. 

UN94 

'  2*120. 

MAXX 

*69. 

CYCMX 

95. 

REZ 

;  2*121. 

MI  NY 

*70. 

CYCPH3 

2*96. 

NODUMP 

i  2*122. 

MAXY 

*71. 

REZFCT 

97. 

UN97 

!  2*123. 

lEXTX 

2*72. 

NTRACR 

98. 

UN98 

;  2*124. 

JEXTY 

2*73. 

NMXCLS 

99. 

UN99 

j  125. 

UNI  2  5 

74. 

BBOUND 

100. 

EVAPM 

!  126. 

UN126 

75. 

UN7  5 

101. 

EVAPEN 

;  127. 

SSI 

76. 

ECK 

102. 

EVAPMU 

i  128. 

SS2 

77. 

NECYCL 

103. 

EVAPMV 

;  129. 

UMIN 

2*78. 

NTPMX 

104. 

EZRH2 

:  *130. 

SS4 

*79. 

GLUED 

105. 

UN105 

131. 

PRTIME 

80. 

UVMAX 

106. 

UNI  06 

!  132. 

EOR 

00 

NTCC 

107. 

UN107 

133. 

EOT 

82. 

UN82 

108. 

UN108. 

134. 

EOB 

2*83. 

IVARDX  • 

109. 

UNI  09 

135. 

EMOR  . 

84. 

T 

*110. 

ROEPS 

*136. 

DXF 

85. 

EMIN 

111. 

UNI  11 

*137. 

DYF- 

*86. 

PMIN 

112. 

UNI  12 

138. 

1 

UNI  3  8 

2*87. 

INTER 

*113. 

FINAL 

i  *139. 

STAB 

88. 

I 

114. 

UN114 

140. 

XIENRG 

89. 

J 

115. 

MBBB 

i  141. 

XKENRG 

90. 

K 

116. 

MBB 

;  142. 

XTENRG 

91. 

M 

117. 

UN117 

1  143. 

UN145 

92. 

N 

2*118. 

NADD 

1  *144. 

DTMIN 

I 

i 
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145. 

UN145 

146. 

EMOT 

2*147. 

JCENTR 

*148. 

RADIUS 

*149. 

BEAR 

1*150. 

EMOB  = 

0  [Last  card  of  SETUP  input  deck.  Not 
included  when  restarting  from  tape.] 

1*R151. 

PK(1) 

should  be  the  same  as  PROB,  problem  number 

*R152. 

PK(2), 

the  restart  cycle. 

*R153. 

PKC3)  : 

when  =  -1.,  program  will  restart  from  tape 

and  do  a  "long"  EDIT  prints  o£  the  pickup 
cycle.  When  =  -2.,  program  will  restart  from 
tape  and  do  a  "short"  EDIT  print  of  the  pick¬ 
up  cycle.  When  =  -3.,  program  expects  to 
restart  from  a  HLPCLM  tape.  Do  not  include 
in  setup  input  deck. 


- 

User- supplied  input-values. 


^  Must  have  a  "1" 

in 

column 

1. 

^  Must  have  a  "2" 

in 

column 

1. 

Must  be  included 

in 

a  restart 

input  deck. 
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APPENDIX  B 

VARIABLES  USED  IN  TRANSPORT  ACROSS  CELL  BOUNDARIES 


Mass 


Energy 


SAMPY^ 

AMPY 


SAMMP. 

1 

AMMP 


SDELETj^ 

DELET 


SDELER. 

DELER 


Radial  Momentum 


Axial  Momentum 


AMUT  AMVT 


AMVR 
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APPENDIX  C 

FLAGS  AND  CONVENTIONS  GOVERNING  INTERFACE  CELLS 


> 

Statements 

Meaning 

RHO(l,M)  =  -1.0 

The  M  location  in  mixed  cell  arrays 
(XMASS,  RHO,  SIE,  FRACTP,  FRACRT, 
SAMMP,  SAMPY)  is  not  in  use. 

NVOID  =  NMAT  +1 

NVOID,  the  number  o£  the  "void” 
package,  is  one  more  than  the  number 
o£  material  packages. 

M  =  MFLAG(K)-100 

RHO (NVOID, M)  =1.0 

Cell  K  contains  the  £ree  sur£ace 
inter£ace . 

MFLAG(K)  =  0 

Cell  K  is  empty  and  is  not  cut 
by  an  inter£ace. 

0  <  MFLAGCK)  <  100 

Cell  K  is  nonempty  and  pure, 
i.e.,  does  not  contain  an  inter£ace. 

MFLAG(K)  =  2 

Cell  K  is  completely  inside  package 

2  boundary  and  contains  only  material 
o£  package  2. 

MFLAGCK)  >  100 

Cell  K  is  mixed,  i.e.,  contains  at 
least  one  inter£ace.  M  =  FLAGCK)-100 
gives  location  o£  quantities  in  mixed 
cell  arrays  £or  cell  K. 

MFLAGCK)  <  0 

Cell  K  was  mixed  on  previous  cycle 
but  is  no  longer  cut  by  an  inter£ace 
and  will  be  re£lagged  at  the  end  of 
PH2. 

MAT(2)  =  3 

The  material  code  number  o£  package  2 
is  3.  The  list  in  EQST  indicates 
that  code  number  3  corresponds  to  iron 

*  1 

♦ 

MATCl)  =  20 

The  material  of  package  1  is  an  ideal 
gas  which  is  given  special  treatment 
in  the  sound  speed  calculation  and  in 
the  strength  phase  of  the  code. 

M  =  MFLAGCK) -100 
RHOCN,M)  =  0. 

The  interface  of  package  N  does  not 
cut  cell  K.  Material  of  package  N 
will  not  be  transported  into  cell  K. 
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Statements 


Meaning 


M  =  MFLAGCK)-100 
RHO(N,M)  >  0. 
XMASSCNjM)  =  0. 

M  =  MFLAGCK)-100 
RHO(N,M)  =  0. 
XMASS(N,M)  >  0. 


Package  N  interface  cuts  cell  K, 
but  cell  K  does  not  yet  contain 
any  material  of  package  N. 

Package  N  interface  has  left  cell 
K,  and  the  material  of  package  N 
[XMASS(N,M)]  should  be  completely 
evacuated  on  this  cycle. 
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